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With  its  low  transmit  power  and  considerable  orbit  distance  from  receivers,  the  global 
positioning  system  (GPS)  waveforms  are  known  to  be  susceptible  to  jamming.  Efforts  to 
date  have  investigated  various  receiver  designs  and  adaptive  antenna  arrays  to  improve  the 
anti-jam  capabilities  of  GPS  position  location,  but  until  this  work  no  extension  of  anti-jam 
capabilities  to  attitude  determination  (determination  of  the  orientation  of  a  body  in  space) 
existed.  This  dissertation  provides  a  comprehensive  investigation  into  the  new  field  of 
robust,  jam-resistant  attitude  determination  using  the  Global  Positioning  System. 

We  present  a  generic  adaptive  antenna  array  and  receiver  design  that  provide  the 
necessary  data  for  both  position  location  and  attitude  determination  in  a  jammed  envi- 
ronment. Using  the  data  this  receiver  provides  we  develop  a  new  maximum  likelihood 
attitude  estimator  (the  MLAE)  that  provides  attitude  determination  capability  even  in 
jammed  environments.  In  order  to  accomplish  this,  new  ways  of  viewing  the  standard 
array  processing  and  direction  finding  tenants  are  presented.  This  leads  to  a  new  in- 
terpretation and  parameterization  of  the  array  response  vector  (spatial  steering  vector) 
and  secondary  data  covariance  matrix.  The  MLAE  optimally  includes  information  from 
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all  satellites  in  view,  and  its  performance  is  shown  herein  to  asymptotically  achieve  the 
Cramer- Rao  Bound  (CRB),  i.e.  the  MLAE  asymptotically  achieves  the  performance  limit 
for  unbiased  estimators. 

As  an  approximation  to  this  estimator,  a  second  estimator  is  developed  that,  at  the 
expense  of  performance,  may  provide  increased  computational  capability.  Two  versions 
of  this  second  estimator  are  presented  that  combine  direction  finding  with  attitude 
estimation. 

A  method  of  optimally  incorporating  data  from  both  GPS  frequencies  into  the 
estimation  is  developed.  This  approach  maintains  the  jammer  mitigation  capabilities  of 
the  previous  estimators  and  allows  for  a  larger  array  spacing,  therefore  increasing  the 
accuracy  of  the  attitude  estimates. 

Simulation  based  performance  results  of  the  various  estimators  are  presented  for 
various  antenna  topologies  and  interference  scenarios.  The  simulation  results  indicate 
that  the  new  algorithms  provide  significant  improvement  over  conventional  attitude 
determination  methods.  In  addition,  the  MLAE  is  shown  to  provide  better  performance 
than  conventional  methods  of  attitude  determination  even  in  unjammed  environments. 
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CHAPTER  1 
INTRODUCTION 

1.1     Motivation 

It  is  well  known  that  adaptive  antenna  arrays  (often  referred  to  as  "Smart  Antennas" 
for  communications  systems)  provide  significant  resistance  to  unintentional  interference 
and  intentional  jamming  for  both  signal  extraction  and  direction  finding.  Global  Position- 
ing System  (GPS)  based  attitude  determination  (AD)  systems  utilize  multiple  sensors  as 
well  to  extract  attitude  through  carrier  phase  differences  between  the  sensors  [1].  However, 
present  AD  algorithms  typically  only  offer  the  jamming  resistance  inherent  in  the  receiver 
(i.e.  gain  from  the  spread  spectrum  waveform  and  perhaps  coupling  with  some  form  of 
INS),  which  for  even  low  powered  jammers  may  not  be  enough  to  prevent  corruption  of 
the  attitude  estimates.  Therefore,  it  seems  appealing  to  exploit  the  similarities  between 
the  two  fields  of  adaptive  array  processing  and  GPS  based  attitude  determination  to  in- 
crease the  anti-jam  capabilities  of  GPS  attitude  systems.  Indeed,  this  work  centers  around 
a  substantially  different  method  of  approaching  the  attitude  determination  task  than  is 
typically  employed. 

The  scenario  that  originally  motivated  this  research  involved  a  medium  to  long 
range  platform  that  must  fly  for  extended  periods  in  a  severe  jamming  environment.  The 
platform  incorporates  an  adaptive  antenna  array  to  provide  GPS  jamming  resistance 
for  position  location  and  navigation.  In  response  to  this  scenario,  a  method  has  been 
developed  that  makes  use  of  an  anti-jam  antenna  array  to  provide  GPS  based  attitude 
information,  as  well  as  the  typical  position  location  data,  even  during  periods  of  jamming. 

Three  distinct  factors  motivated  this  research  into  anti-jam  attitude  determination. 

1.   Threat  from  GPS  jammers.  The  susceptibility  to  jamming  and  inadvertent 
interference  is  high  due  to  the  GPS  satellites'  low  transmit  power  and  considerable  orbital 
altitude.  Even  a  low  power  jammer  several  miles  from  the  GPS  user  can  degrade  accuracy 


or  affect  acquisition.  With  GPS  jammers  being  marketed  as  commercial  items  [2],  one 
could  control  access  to  GPS  over  an  entire  region  by  deploying  several  of  these  low 
cost/low  power  jammers  throughout  an  area  and  activating  them  only  when  desired  to 
deny  the  use  of  GPS. 

2.  Related  research.  Previous  research  in  direction  finding  algorithms  using  antenna 
arrays  for  radar  and  communications  systems  provides  a  substantial  mathematical  basis 
and  a  likelihood  of  tractability  for  the  approaches  developed  in  this  dissertation.  At  the 
macro  level,  attitude  determination  and  direction  finding  are  similar.  For  example,  they 
both  incorporate  multiple  sensors  and  use  carrier  phase  interferometry  to  develop  their 
estimates.  Since  some  direction  finding  algorithms  provide  significant  performance  in 
external  interference  environments,  it  is  natural  to  to  look  to  this  field  with  the  desire  to 
gain  similar  anti-jam  capabilities  for  attitude  determination. 

The  primary  algorithm  developed  in  this  dissertation  does  not  use  direction  finding 
directly.  However,  it  is  not  unlike  many  direction  finding  algorithms  in  the  sense  that 
it  defines  steering  vectors  that  map  unknown  parameters  into  a  model  for  observation 
data  from  the  antenna  array,  and  the  unknown  parameters  are  estimated  by  maximizing 
(or  minimizing)  a  function  of  the  steering  vectors  and  received  data.  Suboptimal  algo- 
rithms developed  later  in  the  dissertation  directly  incorporate  direction  finding,  or  more 
specifically  assume  that  the  directions  are  measured  by  some  method. 

3.  Hardware  trends.  State  of  the  art  systems  are  incorporating  adaptive  antenna 
arrays  with  the  GPS  receivers  to  provide  additional  anti-jamming  resistance  for  position 
location,  providing  much  of  the  hardware  and  design  philosophy  for  attitude  location 
as  well  (see,  for  example,  Falcone  et  al.  [3]).  A  possible  direct  application  of  this  work 
would  be  to  incorporate  the  algorithms  developed  here  into  a  GPS  system  that  already 
incorporates  multi-sensors  for  anti-jam  position  location  and  navigation  in  order  to 
additionally  provide  a  GPS  based  attitude  estimate. 

Responding  to  these  three  motivations  (i.e.  need,  related  research,  and  hardware 
technology  trends),  it  is  natural  to  investigate  a  system  that  makes  maximum  use  of  the 
presence  of  the  anti-jam  sensor  array  by  providing  attitude  as  well  as  position  location. 


1.2     Contributions  of  This  Dissertation 
This  work  draws  on  several  separate  fields  of  research.  The  following  background 
chapter  provides  an  overview  of  the  Global  Positioning  System,  attitude  and  associated 
"conventional"  (i.e.  not  robust  to  interference)  methods  of  attitude  determination  and 
direction  finding  because  all  of  these  fields  were  used  to  develop  new  interference  resistant 
attitude  estimation  methods.  These  results  are  novel  because  external  interference 
mitigation  for  GPS  based  attitude  estimation  is  a  completely  undeveloped  field  [4], 
and  the  results  derived  here  provide  significant  improvement  over  conventional  attitude 
determination  algorithms  in  both  jammed  and  unjammed  environments. 

The  constraints  and  challenges  of  this  task  are  different  from  the  "typical"  antenna 
array  direction  finding  and  communications  applications.  '    To  address  these  challenges 
the  direction  finding  task  is  recast  as  an  "attitude  finding"  problem  by  redefining  such 
standard  and  familiar  terms  as  the  antenna  array  response  vector  (often  called  the 
spatial  steering  vector)  and  data  matrix.  This  provides  increased  insight  into  methods  to 
address  the  problem  and  facilitates  development  of  a  new  maximum  likelihood  attitude 
estimation  method.  This  estimator  is  shown  to  be  asymptotically  unbiased,  consistent, 
and  efficient  (i.e.  asymptotically  achieves  the  Cramer- Rao  Bound).  Another  contribution 
is  the  extension  of  maximum  likelihood  direction/attitude  estimation  to  incorporate  both 
GPS  carrier  frequencies,  as  discussed  in  Chapter  7.  Using  the  method  developed  in  this 
dissertation,  incorporation  of  the  additional  GPS  frequency  provides  a  signal  to  noise  ratio 
(SNR)  gain,  as  opposed  to  methods  in  use  today  in  which  the  SNR  is  reduced  when  the 
additional  frequency  is  used.  Finally,  the  simulation  based  performance  results  provide  a 
quantification  of  the  performance  of  this  work. 


1  Of  course,  caution  is  used  with  the  generalization  of  "typical,"  since  every  application 
has  some  unique  attributes. 


1.3     Dissertation  Overview 
The  remainder  of  this  dissertation  is  organized  as  follows.  Chapter  2  contains  an 
overview  of  the  multiple  research  fields  from  which  this  dissertation  draws.  Chapter  3 
describes  a  generic  receiver  that  provides  the  data  necessary  to  estimate  both  position 
and  attitude,  and  the  differences  between  this  receiver  and  those  typically  employed  for 
either  anti-jam  position  location  or  attitude  determination.  Chapter  4  derives  a  novel 
algorithm  for  attitude  estimation  in  an  external  interference  environment.  This  estimator 
takes  the  form  of  a  maximum-likelihood  estimator.  Chapter  5  presents  several  statistical 
properties  of  this  estimator.  In  Chapter  6,  two  variants  of  suboptimal  algorithms  are 
developed  for  attitude  estimation.  These  algorithms  have  attractive  implementation  and 
computational  features,  but  their  performance  is  not  as  good  as  the  estimator  of  Chapter 
4.  Chapter  7  derives  a  method  of  improving  the  attitude  estimator  of  Chapter  4  by  using 
both  GPS  frequencies.  Chapter  8  presents  quantification  of  the  performance  of  these 
methods  through  simulation  results.  Finally,  Chapter  9  contains  a  summary  of  this  work 
and  a  discussion  of  areas  of  research  related  to  this  work,  but  beyond  the  scope  of  this 
dissertation. 


CHAPTER  2 
BACKGROUND 

2.1     Introduction 

This  chapter  presents  an  overview  of  the  major  fields  which  provide  the  foundation 
for  the  contributions  of  this  dissertation.  Interference  mitigation  for  GPS  based  attitude 
determination  involves  the  synthesis  of  several  different  fields,  so  a  brief  introduction  and 
summary  of  them  is  provided  here. 

This  chapter  is  organized  as  follows.  Section  2.2  provides  an  overview  of  the  Global 
Positioning  System.  Section  2.3  discusses  methods  and  signal  models  used  in  attitude 
determination.  These  two  sections  provide  the  framework  for  the  term  "conventional" 
GPS  based  attitude  determination,  which  is  used  in  this  work  to  describe  the  methods  of 
attitude  estimation  used  today.  Finally,  a  survey  of  direction  finding  is  provided  in  section 
2.4,  with  an  emphasis  on  maximum  likelihood  methods. 

2.2     GPS  Overview 

From  personal,  hand-held  civilian  units  to  integrated  military  systems,  the  Global 
Positioning  System  (GPS)  is  an  important  link  in  the  worldwide  navigation  infrastructure. 
The  Global  Positioning  System  is  designed  to  provide  a  user  with  the  ability  to  calculate 
his  position  in  three  dimensional  space  (i.e.  latitude,  longitude,  and  altitude)  through 
reception  of  signals  from  a  constellation  of  24  dedicated  satellites  orbiting  the  earth. 

The  Global  Positioning  System  conceptually  consists  of  three  segments.  The  space 
segment  involves  the  satellites  in  orbit  about  the  earth.  These  satellites  provide  a  low- 
power  output  signal  that  may  be  received  on  or  above  the  earth.  The  control  segment 
tracks  and  communicates  with  the  satellites  in  order  to  update  their  navigation  informa- 
tion. Finally,  the  user  segment  represents  the  receiver  and  associated  equipment  required 
to  decode  and  interpret  the  signals  sent  by  the  satellites  and  convert  this  into  the  user's 


position  and  velocity.  With  an  impressive  market  penetration,  the  number  of  "users" 
occupying  the  user  segment  is  steadily  increasing. 

The  ability  to  calculate  position  arises  from  measuring  the  range  to  several  satellites. 
If  one  knows  the  range  to  a  fixed  object,  the  locus  of  possible  user  positions  is  on  the 
surface  of  a  sphere  with  its  center  at  the  object  position  and  its  radius  equal  to  the  range 
to  the  object.  Therefore,  if  the  user  knows  his  range  to  three  satellites,  the  positions  of 
the  satellites  (calculated  using  the  satellite's  ephemeris  parameters),  and  has  a  precise 
measurement  of  time,  then  he  can  solve  for  his  position  as  the  intersection  of  the  three 
spheres.  In  practice,  less  precise  clocks  are  used  in  the  user  segment,  requiring  a  minimum 
of  four  satellites  to  obtain  a  position.  Obviously,  errors  and  uncertainty  in  the  ability  to 
measure  the  range  to  the  satellite  directly  affect  the  accuracy  of  the  position  measurement. 

The  GPS  satellites  incorporate  Direct  Sequence-Spread  Spectrum  (DS-SS)  modu- 
lation. This  allows  all  satellites  to  occupy  the  same  spectrum  and,  by  monitoring  the 
code  tracking  loop,  enables  users  to  accurately  calculate  range  to  the  satellite.  Two  levels 
of  accuracy  are  currently  available  through  the  use  of  two  different  spreading  sequence 
architectures.  The  Standard  Positioning  Service,  which  is  available  to  all  users,  incorpo- 
rates a  1023  length  sequence  with  a  chip  rate  of  1.023  MHz  called  the  Coarse  Acquisition 
Code  (i.e.  C/A  Code).  The  Precise  Positioning  Service  (PPS)  uses  the  P(Y)  code  that 
operates  at  a  chip  rate  of  10.23  MHz.  The  P(Y)  codes  have  very  long  natural  periods 
(over  6  months);  however,  they  are  truncated  and  restarted  every  week.  While  the  C/A 
Code  waveform  is  publicly  known  and  easily  obtainable,  the  P(Y)  code  sequence  is  only 
available  to  authorized  users.  Since  knowledge  of  the  spreading  sequence  is  required 
for  demodulation  and  data  detection,  the  PPS  is  essentiaUy  available  only  to  military 
applications. 

The  C/A  waveform  (code  and  data  sequence)  modulates  a  carrier  frequency  of 
1.57542  GHz  (often  referred  to  as  the  Ll  frequency).  Since  several  satellites  are  simultane- 
ously transmitting  at  the  same  frequency,  Code  Division  Multiple  Access  (CDMA)  is  used 
to  allow  unambiguous  satellite  reception.  Each  satellite  uses  a  different  sequence  from  a 
set  of  10th  order  Gold  codes.  Of  the  1025  possible  Gold  codes  available  in  this  family,  37 
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have  been  identified  for  use  on  satellites.  The  P(Y)  code  is  broadcast  on  two  frequencies, 
the  Ll  frequency  defined  above  (in  phase  quadrature  with  the  C/A  code)  and  the  L2 
frequency  of  1.2276  GHz.  At  the  time  of  this  work,  an  effort  is  underway  to  enhance  GPS 
with  additional  frequencies  for  civilian  use  and  a  different  waveform,  denoted  as  the  M 
code,  for  military  and  other  restricted  use. 

The  navigation  message  data  bits  (which  are  modulated  by  the  spreading  sequences 
described  above)  are  transmitted  at  a  rate  of  50  Hz.  The  GPS  navigation  data  consist 
of  synchronization  information,  satellite  ephemeris  data,  and  other  satellite  parameters. 
These  data  are  designed  to  be  packaged  into  30  second  (i.e.  1500  bit)  time  slices. 

An  important  observation  is  that  unlike  many  other  communication  systems  where 
the  only  information  is  the  demodulated  data  stream,  GPS  uses  the  code  timing  and 
tracking  information  to  develop  an  estimate  of  range  to  the  satellites.  This  range,  coupled 
with  the  other  information  obtained  from  the  GPS  data  stream,  allows  a  formulation  of 
user  position. 

2.3     Conventional  Attitude  Determination  Using  GPS 
2.3.1     Attitude  and  Coordinate  Frames 

Attitude  determination  involves  two  coordinate  systems  and  the  transformation  that 
relates  them.  The  local-level  coordinate  system  represents  an  initial,  at  rest,  or  other 
nominal  orientation  of  the  platform  of  interest.  This  frame  is  often  referenced  in  terms  of 
a  "North-East-Down"  (NED)  or  "East-North-Up"  (ENU)  convention.  The  line  of  sight 
(LOS)  vectors  to  the  satellites  are  assumed  known  in  this  coordinate  system.  The  second 
coordinate  system  is  the  antenna  frame  of  reference,  and  represents  the  orientation  of  the 
GPS  antenna  at  a  particular  point  in  time.  The  antenna  frame  and  the  body  frame1   are 
related  by  a  known  transform;  if  one  is  known,  the  other  can  be  easily  found.  For  this 
reason  the  antenna  frame  is  sometimes  referred  to  in  this  work  as  the  body  frame  and  is 


The  body  frame  is  defined  as  the  frame  of  reference  along  principal  platform  axes, 
such  as  along  the  fuselage  and  wings.  Typically  the  body  frame  is  of  higher  interest  from  a 
guidance  and  navigation  perspective. 
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usually  referenced  using  a  "Front-Right-Below"  (FRB)  convention.  The  attitude  is  then 
defined  as  the  relation  between  the  local-level  and  the  antenna  frames  of  reference  at  any 
point  in  time. 

2.3.2     Parameterization  of  Attitude 


Some  common  parameterizations  of  the  attitude  include  Euler  angles,  quaternions, 
and  direction  cosine  matrices.  Throughout  this  work,  this  attitude  is  referred  to  generi- 
cally  as  q,  to  remain  independent  of  any  particular  parameterization.  2    However,  at  times 
a  specific  parameterization  is  required  either  to  clearly  demonstrate  a  point  or  prevent 
ambiguities.  When  specific  parameterizations  of  this  attitude  are  employed  they  follow  the 
notation  and  description  below. 

q  =  the  attitude  of  the  antenna 

£$  =  the  antenna  roll  Euler  angle 

£e  =  the  antenna  pitch  Euler  angle 

£#  =  the  antenna  yaw  Euler  angle 

Q  —  [lii  92, 93, 94 ]T,  the  attitude  quaternion 

Q  =  the  3x3  direction  cosine  matrix  that  transforms  local  level  to  antenna  frame 

The  Euler  angles  represent  successive  rotations  about  the  three  coordinate  axes. 
The  roll  angle,  £$,  represents  rotation  about  the  local  level  x  axis,  the  pitch  angle  £e 
about  the  y  axis,  and  the  yaw  angle  £*  the  z  axis.  The  order  of  rotation  determines 
the  final  orientation.  In  this  work  we  employ  a  "3-1-2"  rotation,  meaning  the  rotations 
are  performed  first  in  yaw,  then  roll,  and  then  pitch.  A  drawback  with  the  Euler  angle 


2  Notice  that  the  unparameterized  attitude  q  in  this  work  is  not  represented  in  bold 
type,  even  though  at  least  three  parameters  are  required  to  specify  the  attitude.  This  is  to 
emphasize  that  when  this  notation  is  used,  the  attitude  may  be  specified  in  any  manner. 


representation  is  the  possibility  of  encountering  singularities;  however,  in  general  they 
provide  an  intuitive,  clear  interpretation  of  the  rotations. 

Another  method  of  defining  attitude  is  with  the  direction  cosine  matrix  Q.  Each  term 
in  the  direction  cosine  matrix  represents  the  cosine  of  the  angle  between  an  axis  of  the 
local  level  frame  and  the  antenna  frame.  For  example,  Q(2, 2)  is  the  cosine  of  the  angle 
between  the  y  axis  of  the  local  level  frame  and  the  y  axis  of  the  antenna  frame  [5].  The 
direction  cosine  matrix  is  orthonormal. 


Qr  =   Q 

IQI   =   i 


(2.1) 

(2.2) 


This  orthogonal  nature  of  the  direction  cosine  matrix  is  very  useful,  since  an  inverse 
rotation  (say,  transforming  from  the  body  frame  back  to  the  local  level  frame)  can  be 
performed  without  inverting  the  matrix. 

The  final  method  of  specifying  attitude  is  with  a  quaternion.  In  order  to  define  the 
quaternion's  ability  to  represent  a  rotation,  first  consider  that  between  any  two  coordinate 
systems  sharing  a  common  origin,  a  single  axis  exists  in  which  measurements  along  that 
axis  are  the  same  in  both  coordinate  systems.  The  relation  between  the  two  coordinate 
systems  (i.e.  the  attitude)  is  then  defined  by  this  axis  and  an  angle  of  rotation  about  it.  If 
we  take  the  axis  as  V  and  the  angle  of  rotation  as  u,  then 


Vt 


V„ 


V, 


where 


|V|  =  1 


(2.3) 


(2.4) 
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and  using  these  the  quaternion  q  may  be  written  as 


vx 

sin(f) 

cosl 

!) 

q  = 


and 

|q|  =  i 

The  identity  quaternion,  which  performs  no  rotation,  is 

0 
0 

0 

1 


(2.5) 


(2.6) 


qidentity  — 


(2.7) 


For  an  excellent  discussion  and  history  of  attitude  and  its  parameterizations,  the  reader  is 
referred  to  Schleppe  [5]. 

2.3.3     Conventional  Phase  Difference  Model 

In  order  to  describe  the  operation  of  a  conventional  attitude  determination  system, 
we  first  describe  the  signal  model  employed  and  the  assumptions  it  incorporates.  Consider 
the  two  GPS  sensor  scenario  of  Figure  2.1.  After  downconversion,  despreading,  and 
integration  the  signal  received  at  sensor  x0  consists  of  a  contribution  from  the  satellite  and 
a  noise  term. 

xQ  =  gei*0  +  n0  (2.8) 

The  response  for  the  second  sensor  may  be  written  in  a  similar  manner. 


Xl  =  geH+o+s*)  +  ni 


(2.9) 


Notice  that  the  second  sensor  also  receives  the  same  satellite  signal  only  shifted  in  time 
(i.e.  phase).  This  phase  difference  is  related  to  the  projection  of  the  baseline  vector  J* onto 


the  LOS  and  the  wavelength: 


8(f>  =  —  v •  d 

A 
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(2.10) 


The  noise  terms  n0  and  n\  are  assumed  uncorrected  and  moreover  are  assumed  to  be 
small  compared  to  the  received  gain  from  the  satellite  signal  g.  The  attitude  determi- 


Sensor  X, 


Sensor  X 


1 


Figure  2.1:  A  two  sensor  interferometric  system.  The  baseline  vector  vector  d  can  be 
determined  by  the  known  LOS  vector  v  and  the  phase  difference  between  the  signals  mea- 
sured at  the  two  sensors  x0  and  x\. 


nation  problem  is  to  estimate  the  baseline  vector  d.  It  is  clear  to  see  that  by  taking  the 
difference  A</>  between  the  phase  measured  at  x0  and  xu  then  the  following  relationship 
between  the  orientation  of  the  vector  d,  the  known  LOS  vector  to  the  satellite  v,  and  the 
phase  difference  may  be  obtained  as  [6] 


A<j>    =    *(xi)-#(xo) 
w    —  (v»d+  k\\ 
w    6</>  +  k(2ir) 

where  $(•)  implies  taking  the  phase  angle  of  (•),  i.e. 


(2.11) 
(2.12) 
(2.13) 


*(•)  =  -3  In(-) 


(2.14) 
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and  w  implies  "is  approximately  equal  to." 

The  reason  that  equation  (2.13)  holds  is  due  to  the  previous  assumption  that  the 
amplitude  of  the  interference  terms  is  much  smaller  than  the  amplitude  of  the  satellite 
terms,  and  therefore  the  satellite  components  dominate  the  phase  of  the  sum  of  satellite 
signal  and  noise.  The  k(27r)  term  represents  the  integer  portion  of  the  phase  difference  be- 
tween the  two  sensors  and  is  called  the  "integer  ambiguity,"  since  it  cannot  be  measured. 
By  using  at  least  three  sensors  (in  a  non-collinear  arrangement)  and  two  satellites  [7], 
the  complete  attitude  of  the  antenna  may  be  found,  once  the  integer  ambiguity  values 
have  been  found.  Several  methods  (see  for  example  Sutton  [8]  and  its  references)  exist  for 
resolving  the  integer  ambiguities. 

Since  the  attitude  estimation  process  makes  use  of  knowledge  of  the  LOS  vectors 
to  the  satellites  in  the  local  level  frame,  a  question  to  consider  is  the  sensitivity  to  error 
in  the  antenna  position  estimate.  It  turns  out  that  small  errors  (within  typical  GPS 
accuracies)  in  position  estimation  will  not  have  a  substantial  impact  on  the  LOS  angles  for 
essentially  every  ground  and  air  based  application.  This  is  due  to  the  large  distance  to  the 
satellites  of  approximately  23,000km  from  the  earth's  center.  Even  for  very  high  altitudes 
of  lOOkft,  the  angular  error  is  negligible. 

The  extension  of  equation  (2.13)  to  multiple  sensors  and  satellites  provides  a  straight- 
forward and  computationally  simple  method  to  estimate  the  antenna  attitude  when  the 
assumptions  of  the  signal  model  are  met.  However,  when  an  interference  source  such 
as  a  jammer  is  present,  the  accuracy  quickly  degrades.  Figure  2.2  indicates  the  level  of 
degradation  a  jammer  can  cause  to  an  attitude  determination  system.  The  parameter 
"mean  total  error"  is  the  average  of  the  total  errors  from  500  realizations  of  the  attitude 
estimation,  where  the  total  error  is  defined  as  the  angular  rotation  between  the  estimated 
attitude  (using  the  attitude  estimator  defined  above)  and  the  true  attitude.3    Notice 
that  even  with  relatively  high  signal  to  jammer  ratios,  the  performance  is  degraded.  The 


Chapter  8  provides  a  mathematical  definition  of  the  mean  total  error. 
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fundamental  thrust  of  this  dissertation  is  the  development  of  algorithms  that  increase  the 
resistance  of  the  attitude  estimator  to  external  interference. 

Mean  Total  Error 
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Figure  2.2:  Mean  Total  Error  Performance  of  the  conventional  attitude  estimator  vs. 
signal  to  jammer  ratio  (SJR),  shown  as  the  solid  line.  The  dotted  line  represents  the  per- 
formance of  the  attitude  estimator  in  an  unjammed  scenario.  Performance  of  the  phase 
difference  approach  to  attitude  estimation  is  degraded  even  when  the  SJR  is  high. 


2.4    Direction  Finding 

The  direction  finding  problem  is  the  determination  of  the  directions  of  arrival  (DOA) 
of  multiple  sources  impinging  on  an  array  of  sensors.  Since  the  problem  may  be  viewed  as 
estimating  how  radiated  energy  is  distributed  over  space,  it  is  often  referred  to  as  "spatial 
spectral  estimation"  [9]. 

Depending  upon  the  application,  DOA  may  be  considered  as  a  single  parameter  (e.g. 
azimuth  angle)  or  two  parameters  (azimuth  and  elevation  angle).  When  the  application 
only  requires  a  single  direction  parameter,  simple  antenna  topologies  like  the  uniform 
linear  array  (ULA)  may  be  employed.  Several  algorithms  make  use  of  this  topology  to 
provide  computationally  efficient  estimates  of  direction.  However,  since  attitude  estima- 
tion is  intrinsically  a  multidimensional  problem  requiring  full  angular  information  to  the 
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satellite  source,  little  insight  can  be  obtained  from  these  single  angle  parameter  algorithms 
for  the  attitude  determination  application.  Therefore,  we  have  limited  the  discussion 
of  this  section  to  those  algorithms  that  can  provide  both  parameters  of  direction  to  the 
source.  This  is  not  to  say  that  the  single  DOA  parameter  algorithms  are  without  merit, 
just  that  they  are  not  applicable  for  attitude  estimation. 

Through  the  last  several  decades  advances  in  direction  finding  theory  have  been  ex- 
tensive, so  much  so  that  even  a  survey  of  the  major  contributions  is  difficult.  One  method 
of  gaining  insight  into  these  advances  is  to  partition  direction  finding  approaches  into 
the  two  classes  of  maximum  likelihood  and  suboptimal.  As  is  often  the  case,  suboptimal 
algorithms  offer  computational  savings,  but  at  the  cost  of  reduced  performance. 

The  first  algorithms  we  consider  are  the  maximum  likelihood  estimators.  These 
optimal  estimators  may  be  further  partitioned  by  considering  the  parameters  estimated. 
The  earliest  works  in  this  field  were  improvements  to  radar  and  sonar  systems  that 
considered  the  angle  and  complex  amplitude  estimation  problem  for  a  single  source  in 
a  jammed  environment  [10].  These  works  focused  on  maximum  likelihood  estimators 
approximated  by  adaptive  monopulse  ratios  (i.e.  adapted  delta  to  adapted  sum  ratio)  [11, 
12].  Although  four  unknowns  exist  in  the  scenario  addressed  by  these  algorithms  (two 
parameters  of  the  complex  gain  and  two  of  the  angle),  only  a  two  dimensional  search 
(over  the  two  angle  parameters)  of  the  maximum  likelihood  metric  is  required,  since  the 
complex  gain  at  the  solution  point  can  be  explicitly  expressed  in  terms  of  the  angle. 

More  recently,  adaptive  array  technology  is  beginning  to  be  applied  to  communication 
systems  to  aid  in  multi-user  wireless  systems  [13].  In  this  extension  of  the  estimation 
problem,  the  unknown  parameters  estimated  are  not  only  the  directions  and  a  complex 
gain,  but  also  an  entire  signal  waveform  for  multiple  signals  impinging  simultaneously  on 
the  array.  For  N  snapshots  of  data  and  L  sources,  there  are  (2 A  +  2)  x  L  real  parameters 
to  estimate  (the  complex  gain  is  incorporated  into  the  unknown  waveform).  As  with 
the  single  source  problem,  the  waveform  parameters  may  be  expressed  in  terms  of  the 
unknown  angles,  resulting  in  a  search  of  the  likelihood  function  over  the  2L  angular 
parameters. 
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Two  distinct  signal  models  have  been  used  in  the  multiple  source  plus  waveform 
estimation  problem,  with  differing  results.  The  "deterministic"  or  "conditional"  maximum 
likelihood  approach  considers  the  source  waveforms  to  be  unknown  deterministic  signals 
in  an  additive  Gaussian  noise  environment.  The  noise  is  often  assumed  temporally 
uncorrelated,  complex  circularly  symmetric,  and  spatially  white.  The  mechanization  of 
the  maximum  likelihood  metric  for  estimating  the  directions  is  to  fit  the  subspace  spanned 
by  the  matrix  of  array  response  vectors  (also  called  steering  vectors)  to  the  sources  to  the 
measurements  in  a  least  squares  sense  [14,15].  Several  direction  finding  algorithms  are  cast 
as  variations  of  the  weighted  subspace  fitting  algorithm  in  Viberg  and  Ottersten  [15]. 

Some  important  results  have  been  found  concerning  the  conditional  maximum 
likelihood  signal  model  assumption.  First,  the  direction  estimates  obtained  from  this 
estimator  are  consistent  [16].  Second,  because  the  number  of  waveform  parameters 
required  to  estimate  increases  with  each  new  snapshot  of  data,  the  waveform  estimates 
are  not  consistent,  and  the  estimator  is  not  statistically  efficient,  i.e.  does  not  achieve  the 
Cramer  Rao  Bound  (CRB),  for  a  finite  number  of  sensors  [16].  However,  as  the  number 
of  sensors  increases,  or  the  signal  to  noise  ratio  increases,  this  estimator  asymptotically 
approaches  the  CRB  [17]. 

The  "stochastic"  or  "unconditional"  maximum  likelihood  approach  is  derived  from  a 
model  where  both  the  signals  and  the  noise  are  stationary,  temporally  white,  zero-mean 
complex  Gaussian  random  processes  [18].  It  has  been  shown  that  the  unconditional 
maximum  likelihood  method  is  statistically  more  efficient  than  the  conditional  ap- 
proach [17, 19].  A  more  computationally  efficient  implementation  that  uses  a  Cholesky 
decomposition  instead  of  an  eigenvalue  decomposition  is  presented  in  Swindlehurst  [18]. 
The  next  class  to  consider  is  estimation  of  the  directions  when  the  signal  waveforms 
are  known  a  priori  to  the  estimation  process.  This  situation  occurs  in  communications 
systems  when  a  known  special  preamble  is  added  to  a  packet  of  data.  In  addition,  the 
GPS  waveforms  fall  into  this  category  once  code  lock  is  obtained.  As  expected,  significant 
performance  gains  may  be  obtained  when  the  signal  waveforms  are  (to  within  a  complex 
gain)  known  as  compared  to  when  they  must  additionally  be  estimated.  In  Li  [20]  and 
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Li  and  Compton  [21]  the  CRBs  for  several  cases  where  the  signal  waveform  is  known  to 
within  a  complex  gain  (and  the  gain  is  either  known  or  partially  known)  are  derived  and 
compared  to  the  CRB  when  the  signal  waveform  is  completely  unknown. 

Finally,  Li  et  al.  show  that  the  estimation  process  of  the  directions  to  the  multiple 
sources  may  be  performed  independently  when  each  of  the  source  waveforms  are  uncor- 
rected with  each  other  [22,23].  This  is  the  so-called  decoupled  maximum  likelihood,  or 
"DEML"  estimator.  The  DEML  is  shown  to  be  consistent  and  asymptotically  efficient  for 
uncorrelated  waveforms.  In  addition,  the  performance  of  this  estimator  does  not  degrade 
as  either  the  number  of  signals  increases  or  the  angular  separation  between  the  sources 
decreases.  The  concepts  presented  in  these  two  works  provide  the  most  direct  applicability 
of  direction  finding  to  the  attitude  estimation  process  derived  later  in  this  work,  since 
the  GPS  P(Y)  waveforms  are  indeed  uncorrelated  and  known,  and  the  number  of  source 
(i.e.  satellite)  waveforms  impinging  upon  the  antenna  array  is  typically  greater  than  the 
number  of  sensors. 

Several  suboptimal  direction  finding  algorithms  have  been  introduced  in  addition  to 
the  maximum  hkelihood  methods  discussed  above,  and  a  few  of  the  more  popular  algo- 
rithms are  now  briefly  discussed.  The  Multiple  Signal  Classification  (MUSIC)  algorithm 
is  based  on  an  eigenvalue  decomposition  of  the  sample  covariance  matrix  [24].  The  algo- 
rithm partitions  this  eigenvalue  decomposition  into  "signal"  and  "noise"  subspaces  and, 
by  observing  that  steering  vectors  to  the  sources  present  will  he  in  the  signal  subspace, 
estimates  directions  by  minimizing  projections  of  steering  vectors  onto  the  noise  subspace. 
Typically,  this  will  involve  a  two  dimensional  search  across  angle.  However,  for  some  array 
topologies  a  significantly  simpler  rooting  algorithm  named  PRIME  may  be  employed  [25]. 
An  algorithm  that  relaxes  the  requirement  for  a  calibrated  array  (one  where  the  steering 
vector  is  known  for  all  possible  direction  of  arrival  angles)  is  the  Estimation  of  Signal  Pa- 
rameters via  Rotational  Invariance  Techniques  (ESPRIT)  [26].  Although  this  is  typically 
an  azimuth  only  algorithm  and  would  therefore  not  meet  our  discussion  criteria,  it  has 
been  extended  to  two  parameter  angle  estimation  by  Swindlehurst  and  Kailath  [27]. 
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Finally,  for  two  excellent  survey  articles  on  a  variety  of  direction  finding  algorithms, 
the  reader  is  referred  to  Godara  [28]  and  Krim  and  Viberg  [29].  In  addition  to  a  summary 
of  the  algorithms,  Godara  [28]  presents  several  implementation  issues  common  in  array 
processing  such  as  array  calibration  (knowledge  of  the  exact  array  response  vector),  failed 
sensors,  and  errors  in  beamformer  weights. 


CHAPTER  3 
AN  ANTI-JAM  GPS  RECEIVER  FOR  ATTITUDE  DETERMINATION 

3.1     Introduction 
In  this  chapter  we  develop  a  functional  GPS  receiver  architecture  that  provides  data 
for  both  position  location  and  attitude  determination  in  a  jammed  environment.  We  apply 
the  following  notation  when  describing  the  receiver.  A  hardware  channel  is  the  portion 
of  the  receiver  from  the  antenna  through  the  analog  to  digital  converter,  containing 
the  mixers  for  down  conversion,  filters,  amplifiers,  etc.  A  satellite  channel  performs  the 
Direct  Sequence-Spread  Spectrum  (DS-SS)  demodulation,  Doppler  demodulation,  and 
integration  for  one  satellite.  Outputs  from  the  satellite  channel  are  used  to  keep  the 
satellite  waveform  in  code  and  Doppler  track,  as  well  as  to  provide  data  for  position 
location.  As  an  example,  a  simple  single-antenna  GPS  receiver  that  can  simultaneously 
track  5  satellites  would  have  (without  incorporating  any  multiplexing)  one  hardware 
channel  and  5  satellite  channels.  Using  this  notation,  we  review  "conventional"  attitude 
determination  and  anti-jam  receivers,  and  then  contrast  them  with  a  proposed  receiver 
designed  to  implement  the  algorithms  of  Chapters  4  and  6. 

3.2     A  GPS  Receiver  for  Anti-Jam  Position  Location  and  Attitude 
Attitude  determination  from  GPS  is  based  on  measuring  the  satellite  time  (phase) 
differences  of  arrival  between  multiple  sensors.  Using  the  terminology  above,  an  attitude 
determination  system  with  M  antennas  that  tracks  L  satellites  would  have  M  hardware 
channels,  and  (again  ignoring  multiplexing)  L  satellite  channels  per  hardware  channel. 
Since  phase  differences,  not  the  absolute  phase  of  the  carrier  sinusoid,  are  used  to  deter- 
mine attitude,  the  M  hardware  channels  may  be  physically  separate  and  independent 
receivers  using  different  clocks,  local  oscillators,  etc.  The  output  phases  from  each  of 
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the  L  satellite  channels  for  all  M  sensors,  with  knowledge  of  the  sensor-to-sensor  base- 
line vectors,  provide  the  data  necessary  for  attitude  determination  in  a  low  interference 
environment. 

Historically,  anti-jam  antenna  arrays  have  operated  as  pre-processors  to  a  GPS 
receiver  [3, 30-32].  The  function  of  the  anti-jam  antenna  is  to  weight  and  combine  the 
voltages  received  at  the  various  sensors  into  one  time  signal,  which  is  then  passed  to  the 
single  hardware  channel.  With  the  interference  attenuated  by  the  adaptive  array,  the 
tracking  loops  following  the  satellite  channels  are  able  to  keep  the  satellites  in  track  and 
provide  data  for  position  location.  Note  that  since  the  anti-jam  antenna  is  often  designed 
as  an  applique,  the  post-antenna  portion  of  the  GPS  receiver  operates  as  if  a  single  sensor 
had  provided  the  time  signal,  even  though  it  is  actually  the  linear  combination  of  data 
from  several  sensors. 

It  is  clear  that  a  pre-processor  type  anti-jam  antenna,  although  it  uses  multiple 
sensors,  does  not  provide  the  data  required  for  attitude  determination.  The  data  signal 
presented  to  the  hardware  chanel  is  the  weighted  combination  of  the  data  collected  by  the 
multiple  sensors.  Even  if  the  weights  are  known,  determining  the  satellite  phase  differences 
across  sensors  from  the  weighted  combination  of  these  sensors  is  an  underdetermined 
problem.  However,  the  receiver  can  collect  the  data  required  for  both  position  location  and 
attitude  determination  by  increasing  the  number  of  hardware  channels  and  incorporating 
the  weight  application  into  the  receiver. 

Figure  3.1  shows  the  top  level  architecture  for  one  possible  GPS  system  that  provides 
adaptive  array-processing  levels  of  anti-jam  performance  for  both  position  location  and 
attitude  determination.  This  system  begins  with  an  array  of  M  sensors.  The  output  of 
each  sensor  is  downconverted  from  the  GPS  carrier  frequency  to  an  intermediate  frequency 
(IF),  using  In-phase  (I)  and  Quadrature  (Q)  downconversion.  The  same  mixing  signals 
(from  the  same  local  oscillator)  are  used  in  each  of  the  M  hardware  channels.  The  signals 
in  each  channel  pass  through  an  anti-aliasing  matched  filter  and  are  sampled  in  I  and  Q. 
For  ease  of  notation  and  discussion,  we  define  a  Vector  Demodulation  block  as  a  col- 
lection of  M  satellite  channels,  one  per  hardware  channel,  with  each  block  demodulating 
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DopplerTrack 
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Location, 
Navigation, . 


Replicate  Beamformers, 
Tracking  updates  for  each 
Satellite  channel 


Attitude  Estimation 


Figure  3.1:  Receiver  block  diagram.  Demodulated  data  are  provided  to  the  attitude  es- 
timator (shown  in  the  circles)  as  well  as  to  the  beamformers.  The  beamformer  output  is 
used  for  position  location  and  updating  the  code  and  Doppler  tracking  loops. 

the  same  satellite  spreading  sequence.  It  is  useful  to  visualize  this  process  as  a  block 
operation  on  all  M  signals  simultaneously,  since  the  same  satellite  demodulation  sequence 
is  applied  to  all  M  signals  from  the  hardware  channels.  The  outputs  of  this  vector  demod- 
ulator are  M  x  1  data  vectors  for  the  various  code  delays  required  to  track  the  satellite 
signal  in  time,  e.g.  an  early,  punctual,  and  late  delay.  More  than  three  code  delays  may 
be  implemented  to  increase  satellite  acquisition  time  or  to  provide  a  quicker  reacquisition 
time  if  the  satellite  signal  is  temporarily  lost.  For  the  multi-sensor  receiver  to  process  L 
satellites,  it  needs  L  vector  demodulators.  The  outputs  of  the  L  vector  demodulators, 
u;,  are  used  separately  for  two  functions,  updating  the  tracking  loops  to  provide  position 
location  information,  and  as  input  to  the  attitude  determination  algorithm. 

Although  the  demodulator  output  vectors  have  received  the  spread  spectrum  process- 
ing gain,  in  a  jammed  environment  the  signal  to  interference  plus  noise  ratio  (SINR)  for 
each  vector  element  may  still  be  too  small  for  satellite  tracking  until  beamformer  weights 
are  applied.  By  projecting  the  M  x  1  data  vectors  onto  the  adaptive  beamformer  weights, 
the  SINR  is  increased  to  an  acceptable  level  to  calculate  the  discriminants  used  for  code 
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and  Doppler  tracking.  Ignoring  the  effects  of  finite  precision  arithmetic,  the  beamformer 
weights  may  be  applied  before  the  satellite  channel  demodulation  (as  with  the  anti-jam 
array  pre-processor)  or  after  satellite  channel  demodulation  (as  is  shown  here)  with  equiv- 
alent results.  In  addition,  the  method  proposed  here  allows  for  different  weights  for  each 
satellite,  improving  tracking  performance  over  a  single  set  of  weights  for  all  satellites. 

The  second  use  of  the  data  vectors  is  for  attitude  determination.  As  discussed  above, 
attitude  determination  requires  the  data  from  multiple  sensors  to  determine  the  time 
(phase)  differences.  Therefore  it  must  operate  on  data  obtained  prior  to  beamforming, 
even  though  this  may  contain  significant  interference  power,  since  as  mentioned  earlier  the 
beamforming  operation  destroys  the  individual  sensor  phases.  The  algorithms  developed 
in  the  next  section  estimate  attitude  even  with  the  interference  present. 

In  summary,  this  section  developed  the  top-level  architecture  for  a  multi-sensor  GPS 
receiver  that  combines  the  anti-jam  benefit  of  adaptive  beamforming  with  development 
of  data  required  for  anti-jam  attitude  determination  algorithms  of  the  next  section.  This 
system  is  like  conventional  attitude  determination  systems  in  that  it  uses  data  from 
multiple  sensors  and  hardware  channels.  However,  it  differs  in  that  adaptive  weights  are 
applied  after  the  satellite  channel  demodulation,  and  that  the  data  passed  to  attitude 
determination  algorithms  are  the  complex  data  vectors  for  each  satellite  and  sensor,  not 
just  the  phase  information.  This  system  is  similar  to  an  adaptive  array  pre-processor 
architecture,  but  differs  in  that  multiple  hardware  channels  (instead  of  a  single  channel) 
are  employed,  and  the  array  weights  are  applied  after  demodulation. 


CHAPTER  4 
MAXIMUM  LIKELIHOOD  ATTITUDE  ESTIMATION 

4.1     Introduction 

In  this  chapter  we  derive  a  new  method  of  for  GPS  based  attitude  estimation.  This 
method  uses  the  receiver  construct  developed  in  Chapter  3.  Motivated  by  Figure  2.2, 
the  desire  is  to  develop  an  attitude  estimator  that  is  jammer  resistant.  This  development 
begins  by  evaluating  the  mathematics  of  jammer  resistant  direction  finding  algorithms, 
and  using  this  background  re-derives  the  approaches  in  terms  of  the  attitude  estimation 
field. 

This  chapter  is  organized  as  follows.  Section  4.2  provides  a  discussion  of  the  signal 
models  and  notation  used  throughout  this  chapter.  In  section  4.3  the  multi-source 
direction  finding  problem  is  reviewed,  to  introduce  the  multi-source  maximum  likelihood 
method.  Specifically  in  this  section  the  task  at  hand  is  to  develop  direction  estimates 
to  each  satellite  in  view.  Section  4.4  presents  some  new  and  modified  concepts  that 
are  required  for  the  new  attitude  estimator,  which  is  developed  in  section  4.5.  Finally, 
section  4.6  provides  a  discussion  of  this  estimator  and  some  graphical  illustrations  of  its 
performance. 
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4.2     Notation,  Preliminaries,  and  Signal  Model 
We  begin  this  section  by  reviewing  the  matrix  and  vector  notation  used  throughout 
this  work. 


Ar 

A" 

A<g>B 

h 

A>0 
Z  6  Caxb 

diag[7i,...,7L] 

qi*q2 


the  transpose  of  the  matrix  A 

the  conjugate  transpose  of  A 

the  Kronecker  (Tensor)  product  of  A  and  B 

the  j  x  j  identity  matrix 

the  j  x  j  matrix  of  zeros 

A  is  positive  definite 

Z  is  a  complex  matrix  of  size  a  x  b 
7!     0      0 


0 


0 


o     o    7L 
quaternion  multiplication 


TR(A)    =    trace  of  A 


In  addition,  the  following  variables  are  consistently  used  in  the  following  roles. 


L    =    the  number  of  satellites  being  tracked 
M    =    the  number  of  sensors  in  the  array 
N    =    the  number  of  snapshots  of  data  available 


(4.1) 


Consider  the  following  popular  baseband  model  for  the  signals  received  by  the  M 
element  array  of  antennas.  The  total  data,  x(t),  received  by  the  M  elements  are  composed 
of  the  contributions  from  each  of  the  L  satellites,  pt(t),  and  the  interference  n(t): 


x(tn)  =  jpiW  +  n(t„),  *»1,...,JV 


(4.2) 


Z=l 
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where  x(t),  pt(t),  and  n(t)  €  CMxl.  We  now  build  up  pt(t)  by  its  components.  Let 
ji  represent  an  unknown  complex  gain  and  0t  the  direction  corresponding  to  the  /th 
satellite  source.  The  direction  0t,  which  is  always  referenced  in  the  antenna  frame,  is 
completely  specified  by  two  parameters.  To  illustrate  this,  consider  an  antenna  located  at 
the  origin  of  the  antenna  coordinate  system.  The  vector  from  the  antenna  to  any  source 
can  be  described  in  spherical  coordinates  by  a  distance  and  two  angular  parameters, 
say  r,  (pi,  and  ip2.  Then  the  two  angles  <pt  and  y?2  comprise  0,  and  the  LOS  unit  vector 
v  =  [vx  vy  vz)T  in  antenna  coordinates  is  found  from: 

cos(</?i)cos(</?2) 

cos(v?1)sin(</?2)  (4.3) 

sin(y>i) 

The  array  response  vector  (spatial  steering  vector),  a(0j),  is  defined  as  the  response  of  the 
array  to  a  signal  impinging  from  direction  0t,  including  both  the  gain  and  phase  response 
of  the  array.  Using  this  definition  the  array  response  vector  can  be  written  in  terms  of  the 
components  of  v  and  a  series  of  vectors  bh  i  =  1, . . . ,  M  from  the  array  reference  point  to 
the  sensor  locations. 


v  = 


a(0)  = 


eJ  a 


e>x*>A> 


(4.4) 


The  signal  from  the  /th  satellite,  8i(t),  contains  the  sampled  spreading  waveform  di(t), 
and  the  sampled  Doppler  eP+*  referenced  at  the  array  reference  point,  defined  as  the 
origin  of  the  antenna  frame. 


si(tn)  =  d,(in)eJUV"  n  =  1,. . .  ,N  (4.5) 

Lastly,  0i(t)  is  the  low  rate  (50  Hz  for  GPS)  data  sequence  that  modulates  the  fast 
rate  spreading  waveform. 
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Using  these  definitions,  the  vector  representation  of  the  signal  received  from  source  I 
at  the  M  sensors  is 


Pi(tn)  m  7*#&)a(0»)*<(*»),  n=l,...,N  (4.6) 

For  this  analysis,  we  assume  that  /?;(£)  is  constant  over  the  period  of  interest,  and 
include  this  constant  in  the  unknown  gain  f\.  If  the  data  bit  is  not  constant  but  its  values 
are  known,  they  can  be  included  in  the  waveform  sequence  sj((„).  Equation  (4.2)  may  now 
be  written  in  a  more  compact  form  as 


L 

=    A(0)rs(tn)  +  n(*n)        n  =  l,...,N  (4.7) 

where 


e  -  [01,02,... A] 

A(0)    =    [aft),  a(fc)f...,  •(«*)] 

r    =    diag[7i,72,...,7i] 

8(*n)      =      [si(tn),S2{tn),...,SL(tn)}T  (4.8) 

and  the  time  varying  components  of  the  satellite  signals  at  time  tn  are  contained  in  s(tn). 

The  signal  portion  of  the  received  data  from  each  source  may  be  interpreted  as  the 
Kronecker  product  of  the  Mxl  spatial  steering  vector  a{0|)  and  the  lxN  temporal  steering 
vector  yi  for  each  source  (we  use  yj  for  the  temporal  steering  vector  to  minimize  confusion 
with  the  column  vector  s(f),  which  contains  information  from  all  sources  at  time  t). 

yi  =  [si(ti),st(t2),...,Sl(tN)}  (4.9) 

This  product  is  the  so-called  space-time  steering  vector  v(fy,y/)  [33], 

v(ft.yi)  =  [yi(f1)aT(0J),  yi(t2)*T(di), . . . ,  yi(tN)aT{6,)]T  (4.10) 
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which  may  be  written  concisely  as 

v(0,,y,)  =  a(0,)®yi  (4.11) 

Using  this  notation  allows  us  to  write  the  MN  x  1  received  data  vector  in  terms  of  the 
space-time  steering  vectors,  complex  gains,  and  additive  interference. 

L 

x  =  ^7,v(0,,y,)  +  n  (4.12) 

where  x  and  n  €  CMN  x  1  as  below. 

x    =    [xT(t1),xT(t2),...,xT(tN)}T  (4.13) 

n    =    [nT(t1),nr(t2),...,nT(tN)}T  (4.14) 

The  model  for  the  interference  is  left  relatively  unstructured,  allowing  for  thermal 
noise,  multiple  jammers,  and  other  interfering  signals.  We  assume  that  the  interference 
noise  waveform  vector,  n(t),n  6  CMxl,  is  zero  mean,  wide  sense  stationary,  and  circularly 
symmetric  complex  Gaussian  with  covariance  matrix  Ra,  where  the  subscript  V  implies 
the  spatial,  i.e.  sensor  to  sensor  covariance.  The  interference  is  uncorrelated  from  all 
satellite  signals.  Moreover,  we  assume  the  interference  is  temporally  uncorrelated.  Since 
thermal  noise  will  exist  in  the  sensors'  data,  it  is  safe  to  assume  that  Ra(0)  is  positive 
definite.  These  conditions  on  the  interference  may  be  concisely  stated  as 

E[n(t)]    =  0  (4.15) 

E[n(<i)nr(tj)]    =  0M  (4.16) 

E[n(U)nH (tj)}    4  RJfr-tj) 

=  Rs(0)<5t4_t.  (4.17) 

R-(0)    >  0  (4.18) 

These  assumptions  are  common  for  analysis  involving  jamming  or  other  interfering 
signals  [21, 22]  and  are  less  restrictive  than  those  of  many  works  on  direction  finding  for 
communication  systems  [16, 34]  which  require  Rs  to  be  a  (scaled)  identity  matrix,  i.e. 
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R.«  =  cr2I.  A  final  assumption  is  that  the  interference  covariance  is  known.  In  practice, 
it  is  not  known,  but  can  be  estimated.  A  method  of  estimating  the  statistics  of  the 
interference  is  presented  in  Appendix  B. 

4.3     Maximum  Likelihood  Direction  Estimation 
In  general,  the  task  is  to  estimate  6L  parameters  from  the  N  samples  of  the  data 
vector  x(t),  these  being  two  parameters  of  the  direction  61,  the  real  and  imaginary  com- 
ponents of  the  complex  gain  7/,  and  the  appropriate  delay  and  Doppler  of  the  temporal 
steering  vector  for  each  source.  However,  the  problem  may  be  simplified  by  assuming 
that  accurate  delay  and  Doppler  estimates  are  provided  from  the  receiver  code  and  phase 
tracking  loops,  providing  yi  as  the  estimate  of  the  temporal  steering  vector  to  source  I. 
This  leaves  \L  real  parameters  to  estimate  from  the  MN  data  samples. 

This  estimation  begins  by  examining  the  likelihood  function.  In  general,  the  space- 
time  likelihood  function  conditioned  on  the  unknown  parameters  0  and  T  is 

/(X,G'r)  =  7r(^)|Rst|eXP  [-(x-  T(0'r)]//R^[x-  T(0,r)]]  (4.19) 

where 

L 

T(e,r)  =  ^7*v(^y0  (4.20) 

(=1 

The  parameters  ©  and  T  that  maximize  (4.19)  are  defined  to  be  the  maximum  likelihood 
(ML)  estimates  of  0  and  T  [35] .  The  space-time  interference  covariance  matrix  Rst 
contains  the  interference  covariances  between  all  sensors  m  =  1, . . . ,  M  at  all  times 
t  =  1, . . . ,  N,  such  that  the  ijth  M  by  M  block  is  E[n(ti)nH  (tj)]    [361- 

By  employing  the  assumption  made  in  (4.17)  that  the  interference  is  temporally 
white,  the  space-time  covariance  matrix  takes  on  a  block  diagonal  structure,  which 
serves  to  decompose  the  argument  of  the  exponential  into  the  sum  of  the  individual  time 
components. 
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Rat  = 


Rs     0     •••     0 

0  Rs  0  0 
0  •••  •.  0 
0      ••■     0     Rs 

Taking  the  negative  of  the  natural  logarithm  of  (4.19),  and  dropping  the  constant 
terms  produces  the  familiar  log  likelihood  ratio  (LLR). 

N 


(4.21) 


LLR  =  £>(«„)  -  A(G)ry(iri)]wR;1[x(in)  -  A(0)ry(*„)] 


(4.22) 


n=l 


Since  maximizing  (4.19)  is  equivalent  to  minimizing  (4.22),  the  estimates  0  and  f  are 
found  as 


N 


0,f  =  argmin  J>(*n)  -  A(0)ry(i„)]HR;1[x(«„)  -  A(0)ry(*n)] 


e,r 


(4.23) 


n=l 


A  significant  factor  in  the  estimation  of  the  directions  to  the  satellites  (and  later,  the 
attitude)  is  the  knowledge  of  the  satellite  waveform.  Since  the  waveforms  are  known,  the 
received  data  may  be  de-spread  and  integrated,  increasing  the  signal  to  interference  ratio. 
This  is  accomplished  when  the  estimated  temporal  steering  vector  is  used  to  de-spread  the 
signal.  De-modulation  of  the  Doppler  and  spreading  sequence  is  achieved  by  projecting 
the  data  from  each  sensor  onto  the  estimated  temporal  steering  (row)  vector  y/  for  each 
source.  This  process  occurs  in  the  "vector  demodulators"  of  Figure  3.1. 

As  shown  in  figure  3.1,  let  uj  represent  this  normalized  projection  onto  the  M  x  N 
data  matrix  X,  which  is  formed  by  reorganizing  the  MN  x  1  received  data  vector  x  such 
that  the  columns  correspond  to  the  time  snapshots  as  shown  below. 


u/ 


Xy/ 


C,H 


(4.24) 
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where 

X    =    [x(tl),x(t3),...,x(tN)]  (4.25) 

N 

-    £>(*»)#(*»)  (4.26) 

n=l 

Each  of  these  projections  are  normalized  by  the  energy  ej  in  the  Zth  source.  This  is  a 
common  normalization  and  is  chosen  so  that  the  power  in  the  integrated  signal  snapshots 
remains  constant,  while  the  signal  to  interference  plus  noise  ratio  (SINR)  increases  linearly 
with  the  number  of  samples  integrated.  That  is,  if  SINR; (TV)  is  the  SINR  in  the  lih 
satellite  component  of  the  data  vector  u;  for  an  integration  of  N  samples,  then 

SINR,(iV  +  1)       JV  +  1 


(4.27) 


SINR,  (TV)  N 

Unless  otherwise  stated,  the  remainder  of  this  work  assumes  that  the  receiver  maintains 
track  on  all  satellites  used  for  direction  or  attitude  estimation  in  both  time  (code)  and 
frequency  (Doppler).  This  allows  the  substitution  y*  =  y,  for  /  =  1, . . .  ,L.  This 
assumption  is  justified  because  the  adaptive  beamformer,  which  follows  the  satellite 
channel,  mitigates  the  effects  of  jammers. 

We  conclude  this  section  by  developing  concise  expressions  for  the  unknown  parame- 
ters. To  review,  the  problem  is  to  determine,  for  each  of  the  L  satellites,  the  2  components 
of  the  direction  vector  and  the  2  components  (real  and  imaginary)  of  the  complex  gain. 
It  can  be  shown  [22]  that  if  the  signal  waveforms  are  uncorrelated,  as  is  the  case  with 
the  waveforms  employed  by  the  satellites,  that  this  multidimensional  estimation  problem 
decomposes  to  a  series  of  decoupled  estimations  of  individual  source  parameters.  At 
the  stationary  point  corresponding  to  the  (as  yet  unknown)  estimate  of  direction  §i,  the 
complex  gain  is  found  to  be 

li=     „  :    '    s     I  l  =  l,...,L  (4.28) 
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and  the  estimator  for  0;  is 

Equation  (4.29)  provides  the  means  to  independently  estimate  the  directions  of  arrival  of 
each  of  the  satellites  in  view,  and  concludes  the  review  of  direction  finding. 

4.4    New  Concepts 
Direction  finding,  in  loose  terms,  states  that  the  information  provided  by  each  source 
(satellite) ,  through  the  differences  in  time  or  phase  of  the  source  signal  as  received  across 
the  sensors,  is  a  direction.  However,  for  the  GPS  attitude  determination  task  consider  the 
new  concept  that  each  satellite  directly  provides  an  (ambiguous)  estimate  of  the  antenna 
attitude  when  the  known  direction  to  the  satellite  in  the  local-level  frame  is  included  in 
the  estimation  process. 

An  Alternate  Array  Response  Vector 

We  begin  by  re-examining  the  array  response  vector  a.  Previously  this  was  defined  as 
a  function  of  the  two  parameter  valued  6,  which  was  derived  from  the  line  of  sight  vector 
v  (in  the  antenna  frame)  to  the  source,  as  shown  in  equation  (4.4).  Indeed,  for  the  last 
30  years  or  so  of  array  processing  the  argument  of  the  array  response  vector  has  been  the 
angular  direction  to  the  source.  However,  another  equally  valid  method  of  parameterizing 
this  expression  is  in  terms  of  the  line  of  sight  vector  in  the  local  level  frame  and  the 
antenna  attitude.  These  two  LOS  vectors  (using  the  subscripts  B  and  LL  to  denote  body 
frame  and  local  level  frame,  respectively)  are  related  by  the  direction  cosine  matrix  that 
transforms  vectors  in  local  level  frame  to  body  frame.  The  direction  cosines  that  comprise 
the  direction  cosine  matrix  are,  of  course,  defined  by  the  body  attitude,  q. 

uB  =  Q(q)  vLL  (4.30) 
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The  new  array  response  vector  parameterized  on  antenna  attitude  is  then  written  as 


&(VLL,q)= 


eJWQCf^u 


(4.31) 


For  the  GPS  attitude  determination  application,  this  method  of  parameterization  is 
preferred  since  the  line  of  sight  vectors  are  known  in  inertial  frame,  and  the  body  attitude 
is  the  desired  quantity.  It  is  important  to  note  that  when  the  antenna  is  actually  at 
the  orientation  incorporated  in  the  second  method,  the  two  array  response  vectors  are 
identical: 

a(0,)  =  a(n,9)  (4.32) 

This  redefinition  of  the  array  response  vector  is  a  critical  step  in  the  development 
of  the  estimator  of  the  next  section.  Notice  that  the  directions  to  the  sources  in  the 
antenna  frame  no  longer  appear  in  the  array  response  vector,  and  supports  the  earlier 
claim  that  this  algorithm  is  not  strictly  a  direction  finding  algorithm.  However,  this  new 
parameterization  introduces  a  complication  not  found  when  parameterized  by  direction: 
attitude  ambiguity. 

Attitude  Ambiguity 

The  attitude  ambiguity  resulting  from  the  new  definition  of  the  array  response  vector 
is  the  manifestation  of  the  fact  that  attitude  cannot  be  uniquely  resolved  when  using 
information  from  only  a  single  satellite  source.  A  useful  visualization  of  this  ambiguity  is 
obtained  by  imagining  the  LOS  vector  from  the  array  reference  to  the  satellite  source  as  a 
"stick,"  and  this  stick  is  "glued"  to  the  array  of  sensors  at  the  array  reference  point.  With 
this  physical  model  the  family  of  possible  attitudes  is  obtained  by  "twirling  the  stick" 
while  keeping  it  pointed  at  the  satellite  source,  as  shown  in  Figure  4.1. 

Mathematically,  this  ambiguity  can  be  obtained  from  the  quaternion  representation  of 
attitude.  Let  q  represent  the  true  attitude  of  the  antenna.  The  locus  of  possible  attitudes 
can  be  determined  by  the  true  attitude,  the  LOS  to  the  satellite  source,  and  a  scalar 
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variable  u  where  -n  <  lu  <  n.  Using  these  factors,  the  ambiguity  in  attitude  from  the  new 
array  response  vector  definition  may  be  represented  as 

a(q>)  =  a((q(u/)*q),  i/)  (4.33) 

where,  from  equations  (2.5)  and  (4.3), 


q(w)  = 


sin(|)i/ 
cos(f) 


(4.34) 


and  *  represents  quaternion  multiplication. 

Fortunately,  this  ambiguity  can  be  completely  resolved  by  using  two  non-colocated 
sources.  This  is  again  consistent  with  standard  GPS  attitude  determination  theory, 
which  states  that  the  minimum  number  of  satellites  required  for  attitude  determination  is 
two  [7].  Appendix  A  proves  this  fact. 

Since  another  common  phenomenon  in  GPS  based  attitude  determination  is  integer 
ambiguity  ,  it  is  worthwhile  to  take  a  moment  to  clearly  delineate  the  difference  between 
the  integer  and  attitude  ambiguities.  The  integer  ambiguity  arises  when  the  sensors  are 
spaced  farther  apart  than  the  Nyquist  sampling  limit.  When  this  occurs,  the  spatial 
spectrum  is  undersampled,  and  possibly  several  discrete  attitudes  would  produce  the  same 
phase  differences  across  the  antenna  array.  This  is  the  same  phenomena  as  the  antenna 
design  term  grating  lobes.  As  shown  above  and  in  Figure  4.2,  the  attitude  ambiguity  in  the 
array  response  vector  is  actually  a  continuum  of  possible  attitudes,  not  a  finite  discrete  set 
as  in  the  case  of  the  integer  ambiguity. 

4.5     Maximum  Likelihood  Attitude  Estimation 
With  the  re-parameterization  of  the  array  response  vector,  we  now  proceed  with 
the  development  of  the  maximum  likelihood  attitude  estimator.  The  observables  for 
this  algorithm  are  the  vector  outputs  of  the  satellite  channels,  i.e.  the  demodulated 
and  integrated  data  vectors  u;.  There  are  two  reasons  this  algorithm  only  considers  the 
demodulated  and  integrated  data  vectors.  First,  the  spreading  sequence  and  Doppler 
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Array  Reference  '* 

Figure  4.1:  Visualization  of  the  possible  (ambiguous)  attitudes  corresponding  to  a  single 
source.  The  possible  array  attitudes  can  be  found  by  "twirling"  the  LOS  vector,  which  is 
assumed  affixed  to  the  array,  while  keeping  it  pointed  at  the  satellite  source. 

demodulation  is  typically  implemented  in  hardware  [37].  Second,  the  GPS  signal  is  too 
weak  to  evaluate  before  significant  integration. 

Recall  from  equations  (4.12)  and  (4.24)  that  uh  I  =  1, 2, . . .  L  contains  contributions 
from  all  satellites  and  interference.  We  denote  the  signal  contribution  from  the  /th  satellite 
as  7,1  and  the  interference  (including  jammers,  thermal  noise,  and  the  multiple  access 
interference  from  other  satellites)  as  wj,  so  that 


U(  =  Zt  +  W; 


(4.35) 


where 


(4.36) 


and 


W| 


Sl 


[n(t1),n(t2),...,n(tN)]+    £    7pa(i/p  ,q)yf 

p=i,p^( 


yf 


(4.37) 
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Figure  4.2:  Possible  roll  and  pitch  Euler  angles  arising  from  the  ambiguity  in  the  new  ar- 
ray response  vector  (yaw  angles  not  shown).  The  true  roll  is  25  degrees,  and  the  true  pitch 
is  15  degrees,  as  indicated  by  the  diamond.  The  locus  of  possible  attitudes  is  determined 
by  the  LOS  to  the  source  and  the  true  attitude.  A  different  source  LOS  would  have  a 
different  ambiguity. 

and  q  represents  the  actual  (true)  attitude  of  the  antenna  array.  As  was  assumed  pre- 
viously, the  signals  are  assumed  in  track  such  that  yj  is  essentially  yj.  It  is  clear  that 
u;  is  a  consistent  estimator  of  zt.  Just  as  with  the  thermal  noise  and  jammer  signals  we 
consider  w*  to  be  a  zero  mean,  circularly  symmetric,  complex  vector  of  Gaussian  random 
variables.1 

E[w,]    =    0  (4.38) 

E[wjwf]     =    0  (4.39) 

Using  the  these,  we  define  the  space-satellite  data  vector,  U,  of  the  received  data  as 


1  Using  this  model  for  the  multiple  access  interference  is  similar  to  approaches  used  for 
bit  error  analysis  of  multi-user  communication  systems. 
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Ul 

u-    U! 

Z  and  W,  the  components  of  U,  are  similarly  defined  as 


Zl 

z2 


•             -1 

Wi 

w  = 

w2 

wL 

(4.40) 


(4.41) 


(4.42) 


such  that 

U  =  Z  +  W  (4.43) 

Finally,  let  Rss  denote  the  ML  x  ML  space-satellite  covariance  formed  from  W  as 


RS3  =  E  [WW"] 


(4.44) 


As  with  the  direction  finding  example  earlier,  the  attitude  estimation  problem  is  cast 
in  terms  of  the  log  likelihood  function.  However,  now  the  LLR  is  a  parameterized  by  the 
antenna  attitude  and  the  complex  gains,  and  as  with  any  likelihood  function  the  true 
parameter  (i.e.  the  true  attitude  q)  is  replaced  with  a  variable  (q),  the  parameter  to  be 
estimated: 


LLR  =  [U  -  Z(T,q)]"  R"i[U  -  Z(I\<7)] 


(4.45) 


Using  the  results  of  Appendix  B,  the  ML  x  ML  space-satellite  interference  covariance 
matrix  asymptotically  reduces  to  a  block  diagonal  structure  ofMxM  blocks: 
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*V«j 


(4.46) 


Cx      0      •••      0 

0     C2     0       0 

0      •••     '••      0 

0      •••      0     CL 

Therefore  as  N  becomes  large,  the  LLR  asymptotically  approaches  the  following  sum  of 
terms: 

L 

LLR  ~  J}u,  -  jMvt^fCT^Ui  -  7ja(i/,,g)]  (4.47) 

ui 

Notice  that  we  have  reparameterized  the  LLR  in  terms  of  the  known  satellite  LOS 
directions  in  the  inertial  frame,  the  desired  attitude  q,  and  now  the  unknown  complex 
gains  7j. 

As  with  the  direction  finding  application,  a  closed  form  expression  for  7j  at  a  station- 
ary point  may  be  obtained  by  setting  the  partial  derivative  of  (4.47)  with  respect  to  7; 
equal  to  zero  and  solving  for  7;  in  terms  of  q. 

aH(vhq)C^ui 


li 


(4.48) 


Substituting  this  into  equation  (4.47)  produces  in  equation  (4.49)  the  Maximum  Likeli- 
hood Attitude  Estimator  (MLAE): 

L 

q  =  argmin]T[u(  -  7/a(i/(,g)]//C(-1[u(  -  7ta(i/,,g)]  (4.49) 

One  way  to  view  this  expression  is  that  it  produces  a  "metric"  value  at  every  possible 
attitude,  and  the  attitude  with  the  lowest  metric  value  is  chosen  as  the  estimate. 

We  now  present  three  alternate  forms  of  the  MLAE.  The  first  alternate  method  of 
expressing  the  MLAE  is  found  by  expanding  the  summation  above: 


q  =  min  £  ufCf »«,  -  7,*a(i/,,  qfC^m 


1=1 


7(ufC(  la(l>|,  q)  +  7^7j*a(i/,,  q)HCjlaL{vuq)    (4.50) 
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and  substituting  in  the  expression  for  the  estimate  of  the  complex  gain,  7,  from  (4.48)  to 
produce 


L 


A  simpler  form  of  the  expression  may  be  written  using  the  shortened  form  of  the  array 
response  vector: 

ai  =  a(i/,,g)  (4.52) 

A  "whitened"  array  response  vector  may  be  written  in  shortened  form  as  well: 

a,  =  Cr^afa,?)  (4.53) 

where 

C,    =    C;/2C(1/2  (4.54) 

c;/2  =  (&/>)"  (4.55) 

Using  these  definitions,  (4.51)  may  be  rewritten  as 

l/2;w„\i2 

(4.56) 

arg  max  J]  uf  CT^P^ jPf^Ui  (4.57) 


A  |uf  Cr1/2a,(g)|2 
q    =    argmax>  '  ,  >,,      ' 

9     M  W^l 


L 


9       (=1 


where 


P^  -  WWM)  (4-58) 

Finally,  a  very  compact  expression  for  this  estimator  for  q  may  be  written  as 

L 

^argmin^TRpjC,-1]  (4.59) 

9    ;=i 

where 

P,  =  [u,  -  iia(vi,  q)]  [m  -  7,a(i/j,  q)\H  (4.60) 

4.6     Discussion 
The  most  salient  observation  of  this  estimator  is  that  the  likelihood  metric  is  defined 
as  a  function  of  the  antenna  (body)  attitude.  Directions  to  sources  in  the  antenna 
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frame  are  not  a  part  of  this  estimator.  Regardless  of  how  parameterized  (Euler  angles, 
quaternion,  rotation  matrix,  etc.),  the  attitude  that  minimizes  (4.49)  is  the  maximum 
likelihood  estimate  of  the  antenna  attitude,  and  since  the  space-satellite  covariance  matrix 
is  (for  large  sample  numbers)  block  diagonal,  the  estimator  decomposes  to  a  sum  of  terms, 
each  of  which  contains  contributions  from  a  single  satellite. 


Yaw  degrees 


Pitch  degrees 


Figure  4.3:  Contribution  to  the  value  of  the  maximum  likelihood  attitude  estimator  of 
equation  (4.49)  from  the  first  of  seven  satellites  in  view.  This  satellite  provides  little  in- 
formation about  yaw,  and  the  most  in  the  positive  pitch  -  negative  yaw  to  negative  pitch 
positive  yaw  dimension. 


To  provide  insight  into  the  operation  of  the  estimator,  examine  these  contributions  in 
a  neighborhood  near  the  true  attitude.  When  the  LOS  to  an  interference  source  is  close 
to  that  of  the  satellite,  intuitively  the  contribution  of  this  satellite  to  the  final  solution 
should  be  small.  Indeed,  this  is  the  case  with  this  estimator.  In  the  neighborhood  of  the 
true  attitude,  the  contribution  to  the  metric  of  satellites  close  to  an  interference  source 
varies  much  less  than  the  contribution  from  those  with  a  large  spatial  separation  from 
interfering  signals.  Visually,  this  causes  the  metric  to  appear  flatter  across  attitude,  and 
therefore  contribute  less  to  the  final  attitude  estimate.  The  shape  of  the  likelihood  metric 
contribution  from  a  particular  satellite  is  not  uniformly  steep  or  flat  across  all  attitude 
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Figure  4.4:  Contribution  to  the  value  of  the  maximum  likelihood  attitude  estimator  from 
the  fifth  of  seven  satellites  in  view.  This  satellite  provides  little  information  in  the  positive 
pitch  -  negative  yaw  to  negative  pitch  -  positive  yaw  dimension,  incidentally  the  dimension 
satellite  in  which  satellite  one  provided  the  most  information. 


dimensions,  but  varies  with  the  amount  of  attitude  information  that  satellite  provides  in 
each  dimension.  For  example,  in  a  non-jammed  scenario  a  satellite  directly  above  a  level 
antenna  array  provides  little  information  about  the  antenna  yaw  angle.  Similarly,  the 
location  of  jammers  may  decrease  the  amount  and  type  of  attitude  information  a  satellite 
may  provide.  In  effect,  the  estimator  is  weighting  each  satelhte's  contributions  to  the 
attitude  estimate  by  the  amount  and  type  of  information  that  they  are  able  to  provide. 
This  phenomena  may  be  observed  in  Figures  4.3,  4.4,  4.5,  and  4.6.  Figures  4.3,  4.4, 
and  4.5,  depict  the  contribution  of  three  individual  satellites  to  the  likelihood  "metric"  as 
a  function  of  yaw  and  pitch.  Only  the  two  attitude  dimensions  yaw  and  pitch  are  shown 
for  graphical  reasons;  for  the  plots  the  metric  is  evaluated  at  the  true  roll  angle.  The  true 
attitude  in  Euler  angles  is  [0  0  0],  i.e.  the  local  level  frame.  Note  that  the  individual 
satellite  contributions  contain  regions  where  they  provide  high  quality  information  (i.e.  a 
steep  slope  across  a  dimension)  and  regions  which  contribute  little  information  (little  slope 
or  fiat).  When  the  individual  contributions  are  combined,  one  satellite's  weakness  may 
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Figure  4.5:  Contribution  to  the  value  of  the  maximum  likelihood  attitude  estimator  from 
the  seventh  of  seven  satellites  in  view. 

be  compensated  by  another's  strength,  producing  the  steep  slope  in  all  dimensions  easily 
observable  from  Figure  4.6. 

The  nature  of  the  likelihood  metric  has  implementation  benefits  as  well.  This 
separable  by  satellite  structure  provides  a  framework  for  parallel  computations,  where 
for  a  given  attitude  under  test,  the  terms  of  the  likelihood  for  each  satellite  could  be 
computed  simultaneously  on  multiple  processors.  In  addition,  should  the  receiver  lose  lock 
on  a  particular  satellite  signal,  its  contribution  to  the  metric  can  easily  be  removed. 

4.7    Comments  on  Searches 
There  are  many  ways  the  two  dimensional  search  of  (4.29)  or  the  three  dimensional 
search  of  (4.49)  can  be  implemented.  In  general,  the  likelihood  surface  does  not  monotoni- 
cally  converge  to  a  single  global  minimum,  but  may  have  several  local  minima.  However, 
if  a  reasonably  accurate  initial  estimate  of  the  antenna  attitude  is  available,  it  can  be  used 
to  provide  a  starting  point  for  a  limited  search  over  an  uncertainty  region  that  contains 
only  one  minimum.  One  method  of  searching  this  surface  is  to  evaluate  (4.29)  or  (4.49)  at 
several  points  in  a  coarse  grid  spanning  the  uncertainty  region.  The  point  with  the  largest 
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Figure  4.6:  Total  value  of  the  metric  including  contributions  from  all  satellites  in  view. 
Notice  that  the  areas  of  weak  information  from  any  particular  satellite  have  been  filled  in 
by  the  contributions  of  other  satellites. 

likelihood  value  is  then  chosen  as  the  center  of  a  finer  grid  search,  where  the  fine  grid  size 
is  chosen  from  the  required  angular  resolution  for  the  particular  application. 

Another  method  is  to  use  a  variation  of  the  alternating  maximization  proposed  by 
Ziskind  and  Wax  [14]  (for  further  reference  see  Press  et  al.  [38]).  For  a  two  dimensional 
direction  finding  search  the  process  is  as  follows:  First  6X  is  fixed  at  some  initial  estimate 
and  the  LLR  is  maximized  with  respect  to  6y  across  the  uncertainty  area.  Then  9y  is 
fixed  at  the  value  that  maximized  the  likelihood  function  and  the  LLR  is  maximized 
with  respect  to  0X.  This  iterative  method  alternates  movement  throughout  the  two 
dimensional  likelihood  space  (parameterized  by  6X  and  6y)  along  the  Qx  and  Bv  axes,  and 
will  converge  to  a  local  minimum.  The  extension  to  a  three  dimensional  attitude  search  is 
straightforward.  Again,  if  the  uncertainty  area  is  reasonably  small  and  contains  the  true 
global  maximum,  then  this  method  will  converge  to  it.  If  no  a  priori  information  exists  to 
provide  an  initial  estimate  of  attitude  (and  therefore  no  initial  estimate  of  6  for  the  search 
of  equation  (4.29)),  then  the  problem  is  significantly  more  tedious.  In  this  case  use  of  a 
more  exotic  method  of  searching,  such  as  genetic  algorithms  may  be  employed  [39,40]. 


CHAPTER  5 
PROPERTIES  OF  THE  MAXIMUM  LIKELIHOOD  ATTITUDE  ESTIMATOR 

5.1     Introduction 

Chapter  4  derived  a  new  method  of  estimating  antenna  attitude  using  the  demodu- 
lated (despread)  data  vectors  for  the  satellites  in  view,  the  maximum  likelihood  attitude 
estimator  (MLAE).  It  showed  that  by  reparamterizing  the  array  response  vector  in  terms 
of  attitude  and  not  antenna-frame  direction,  the  estimator  decomposed  into  a  series  of 
components,  one  per  satellite  source,  parameterized  over  the  antenna  attitude. 

This  attitude  estimator  is  fundamentally  different  from  the  direction  finding 
works  [22, 23],  where  each  source  direction  estimation  is  decoupled  from  another,  even 
though  in  both  the  attitude  determination  and  direction  finding  cases  the  source  wave- 
forms are  known  and  uncorrected.  However,  in  this  chapter  we  show  that  many  of  the 
desirable  statistical  properties  of  the  direction  finding  estimators  apply  to  the  new  attitude 
estimator  as  well.  Specifically,  in  this  chapter  we  show  that  the  MLAE  is  asymptotically 
consistent,  unbiased,  and  efficient. 

Recall  from  Chapter  4  that  the  MLAE  attitude  estimate  (from  equation  (4.56),  which 
is  repeated  here  for  convenience)  is  the  attitude  that  maximizes  a  sum  of  terms  obtained 
by  despreading  the  sensor  data  matrix  with  each  satellites  spreading  waveform: 

L 
q  =  max£  uf  Cr1/2Pa,(9)Cr1/2u(  (5.1) 

where 


9      Ul 


_  a,(<?)af(g) 

*™  "  WHfiSB  {  ] 

Now  consider  the  right  side  of  equation  (5.1)  as  the  function  G(q)  of  the  unknown 
attitude: 

G(q)  =  f:u»C;^l{q)C^\  (5.3) 
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G(q)  is  simply  the  value  of  the  estimator  "metric,"  parameterized  on  the  various  hypothe- 
sized attitudes  of  the  GPS  antenna  array,  as  discussed  in  Chapter  4. 

This  chapter  is  organized  as  follows.  First,  using  the  function  G(q)  defined  above,  this 
chapter  proves  that  the  MLAE  is  a  consistent  estimator.  Consistency  is  a  requirement  for 
the  development  of  the  remainder  of  the  chapter.  To  facilitate  proof  that  the  MLAE  is 
unbiased  and  efficient,  a  series  of  lemmas  on  the  asymptotic  properties  of  various  aspects 
of  the  MLAE  are  presented.  Using  these  lemmas,  this  chapter  then  presents  two  theorems 
that  show  that  the  MLAE  is  asymptotically  unbiased  and  statistically  efficient.  The 
chapter  concludes  with  a  summary  and  a  brief  discussion  of  the  validity  of  asymptotic 
properties. 

5.2     Consistency 
A  consistent  estimator  has  the  desirable  property  that  the  estimation  error  decreases 
as  the  number  of  data  samples  increases.  For  the  MLAE,  this  corresponds  to  the  number 
of  samples  used  for  integration,  N. 
Theorem  1    The  MLAE  is  a  consistent  estimator  for  antenna  attitude.    That  is, 

J™  q  =  q  (5.4) 

Proof:     Consider  the  contribution  of  only  the  first  satellite  to  G(q)  in  equation  (5.3), 
and  define  this  contribution  as  Gi(q): 

G1(q)  =  n^C;1/2P,l{a)C;1/2ul  (5.5) 

Now  from  Chapter  4  uj  is  a  consistent  estimator  of  7^(1/!,  §)  where  q  represents  the  true 
attitude.  So  asymptotically,  (5.5)  will  converge  to  the  maximum  of 

s,H(vuq)CTlSi(uuq)2 

h 


«tW)cr1«fa,S) 

^{u^C-M^q)  (5-6) 


or  using  the  whitened  array  response  vector  notation,  the  maximum  of 


•2 

&(vi,q)B&(vi,q) 

171 '    ^.jWn.g)  (5'7) 
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where  a.{vuq)  is  defined  in  (4.53)  One  could  look  at  equation  (5.7)  and  ascertain,  by 
the  Cauchy-Schwarz  inequality,  that  the  attitude  that  maximizes  G^  is  q.  However,  this 
would  not  be  entirely  correct.  Indeed,  q  would  be  one  attitude  to  maximize  Gx.  But 
since  the  array  response  vector  is  ambiguous  when  parameterized  in  attitude,  an  entire 
family  of  attitudes  will  maximize  Gx.  Let  qi(u>i)  represent  the  family  of  attitudes  that 
asymptotically  maximize  (5.7).  Then  using  (4.33)  and  (4.34),  q^)  can  be  written  as 

qi(wi)=qi(wi)*q  (5.8) 

where  •  represents  quaternion  multiplication. 

Now  consider  the  remaining  satellites'  contributions.  In  a  similar  fashion  as  with 
satellite  one,  the  family  of  attitudes  q2(uj2)  that  asymptotically  minimize  the  second 
through  Lth  satellite's  contribution  to  (4.49)  can  be  found  as 

qi(wi)=qi(o;i)*q        i  =  2,...,I  (5.9) 

For  the  L  satellites,  there  now  exist  L  families  of  attitudes  that  separately  minimize 
each  satellites  contribution  to  (4.49).  In  general,  the  parameter  value  that  minimizes 
a  sum  of  functions  is  not  necessarily  the  value  that  minimizes  any  particular  function. 
However,  for  the  attitude  determination  case,  Appendix  C  states  that  each  of  these 
families  of  minima  intersect  in  a  single  point,  the  true  attitude  q.  So  asymptotically,  the 
estimated  attitude  is  the  true  attitude,  and  therefore  the  MLAE  is  consistent.  ■ 

5.3     Lemmas  on  the  MLAE 
We  begin  by  considering  the  function  G{q)  of  equation  (5.3)  near  the  estimated  atti- 
tude q.  A  popular  method  for  analyzing  G(q),  and  therefore  the  estimator  performance, 
is  through  a  Taylor  series  expansion  of  the  gradient  vector  g'(q)  near  q  (see,  for  exam- 
ple [16]).  The  gradient  vector  is  a  3  x  1  vector,  since  there  are  three  independent  attitude 
parameters  (for  example,  the  three  Euler  angles  fc,  fe,  and  £*).  From  (5.1),  q  maximizes 
G(q),  and  therefore  the  gradient  g'(q)  at  this  stationary  point  is  zero.  If  the  difference 
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between  the  estimated  and  true  attitude  is  small,1   then  the  higher  order  terms  in  the 
Taylor  series  expansion  may  be  ignored,  leaving: 


0  =  g'(q)=g'(q)  +  g"(q)(q-°q)  + 


(5.10) 


such  that 


where 


(q  ~  °q)  =  ~    E"(q)       g'(?) 


-l 


S'(q)  = 


-lT 


9G(q)        dG(q)       dG(q) 
dqi  dq2  dq3 


and 


g"(q)  = 


( d2G(q)  d2G(q)  d2G(q)\ 

dqidqi  dqidq2  9qidq3 

92G(g)  d2G(q)  d2G(q) 

dq2dq\  9q2dq2  9q2dq3 

92G(q)  d2G(q)  d2G(q) 

\dq3dq2  dq3dq2  dq3dq3  J 


(5.11) 


(5.12) 


(5.13) 


and  <7i,  q2,  and  q$  are  the  three  attitude  parameters.  We  now  develop  four  lemmas 
describing  the  terms  of  equation  (5.11). 

Lemma  1    The  gradient  vector  g'(q)  of  equation  (5.12),  evaluated  at  q  =  q,  asymptotically 
may  be  written  as 

g'(q)c±2Re{  (5.14) 


EtfD"(0Pj(5)w«} 


Proof:     This  lemma  is  proven  by  straightforward  evaluation  of  the  terms  in  (5.12). 
We  begin  by  evaluating  the  partial  derivative  of  G{q)  with  respect  to  the  ith  attitude 
parameter. 


9G(q) 
dqi 


=  2  Re  K>f  cr^'p^^a^/^^cr1^ 


i=i 


9=9 


where 


iKq)  =  (a?(q)*i(q)r1a?(q) 


(5.15) 


(5.16) 


1  The  error  is  small  if  the  integration  time  (i.e.  ./V)  is  sufficiently  large  since  the  MLAE 
is  a  consistent  estimator,  as  previously  shown. 
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and  we  define 

4(1)  =  -^-  (5.17) 

and 

dt(l)  =  Cr1/ad«(Z)  (5.18) 

and  have  used  the  fact  that  the  derivative  of  the  projection  matrix  is  [41] 

^Pi,(«)  =  PiMMl)3l(q)  +  (Pi(,)di(0Sj(9))  (5-19) 

In  (5.12)  and  the  remainder  of  this  chapter  we  use  the  shortened  array  response  vector 
notation  of  (4.52).  Using  the  method  of  [22,23],  the  second  u;  in  (5.15)  may  be  replaced 
by  its  asymptotic  value  jian(q). 


2Re|x:u^Cr1/2Pafi(?)di(/)a(t(g)Cr1/%a,(g)|  (5.20) 


gggB 

dq{ 

=    2Re|^Ufcr1/2Pafi(?)di(/)a/(^)7ia((9)|  (5.21) 

=    2Re|^u^Cr1/2Pafi(?)di(/)7i|  (5.22) 

and  observing  that  Pf  „  C~ll2a^{q)  is  zero,  equation  (5.22)  may  be  equivalently  written  as 


■id) 

L 


^»     .    2Re{^(u,-7la,(§))"c,-"¥i(„)a(/)7,J  (5.23) 

=    2Itej£w*C,-''2Pf;(.)a,(/)7,l  (5.24) 

where  Wj,  defined  in  (4.37),  may  be  written  as 

w,  =  uj-7,aj(g)  (5.25) 

and 

-,-1/2 


w,  =  C,   '  wj  (5.26) 


Therefore,  the  zth  term  of  the  gradient  vector  may  be  asymptotically  written  as 

dG(g) 
dqi 

Now  define  the  m  x  3  matrices 


2  Re 


E^f(0Pf(?)wj 
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(5.27) 


D(Z)  =  [d1(0d2(/)d3(0] 


and 


D(i)=  [di(0d2(0d3(0 
Using  (5.27)  and  (5.29),  the  entire  gradient  vector  g'(q)  is  found  as 


g'(?)    ± 


dG(q)       dG(q)       dG(q) 
dqi  dq2  dq3 


-  2Re{t^H(lKc^] 


(5.28) 


(5.29) 


(5.30) 


which  concludes  the  proof. 

Lemma  2   The  Hessian  matrix  g"  of  equation  (5.13)  asymptotically  may  be  written  as 


g"(g)  =  -2Re|x;|7i|2D^(/)P^D(o| 


(5.31) 


Proof:     This  lemma  is  proven  by  evaluation  of  the  terms  in  (5.13).  We  begin  by 
considering  the  partial  derivative  of  (5.15)  with  respect  to  the  itth  attitude  parameter. 


d2G(g) 
dqidqk 


=  2  Re   5>?c; 


1/2 


1=1 


dqk 


P£d,-(0aJl  c,-1/2u, 


(5.32) 


where 


Now 


examine 


dlk 


a/  =  (af  5j)  a/*  (5.33) 

P±d\(/)aJJ ,  the  term  to  be  differentiated.  Prom  the  definition  of  a 


P^ ,  this  may  be  rewritten  as 
8 


dq; 


-[PiMitf]   =  g- [hcosJ  -  p^osJ] 


d_ 

dqk 
d 

dqk 


[*-P«,*] 


(5.34) 
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where  we  define 


st 


*  =  d,-(/)aJ 


(5.35) 


Using  the  chain  rule  to  evaluate  the  derivatives  of  the  second  term  in  (5.34),  we  get 

=     (4(1))^  if +  3,(0(1^  (5.36) 


dqk 


where  ()'k  implies  differentiation  with  respect  to  the  kth  attitude  parameter.  Using  (5.19) 
and  (5.36),  the  second  term  of  (5.34)  is  found  as 

C|  ft  r\ 

dqk  dqk  0ft 


=    P£,dfc(/)a|di(Z)a;+(ait)    df  (0P£d,-(0a| 
+PS|(d<(0)^alt  +  Pi,d1(0(ift)|k 

The  derivative  of  the  psuedo- inverse,  aj,  is  found  by  direct  evaluation  to  be 
=    (a* a,)"1  af  (I)  -  (if  a,)"2  (df  (0*  +  if d*(l))  if 

=  (if^-^dfcopi-^iaji/) 


(5.37) 


(5.38) 


Substituting  (5.38),  (5.37),  and  (5.36)  into  (5.34)  produces  the  desired  partial  derivative: 

-  (^)*V(0P4*(0<+ W*)"xPta,(iA(iy"P4 

+  (a//ai)"1Pa-,di(/)a//d,(/)a/t    (5.39) 


which  may  be  consolidated  and  rewritten  more  simply  as 

a 


% 


P^Oaj]  =  -  (i})*df  (0Pi*(0iJ  +  P£*a 


(5.40) 
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by  defining  #2  to  be 


*2  =  (di(O)^a!  -  afc(I>|}di(04  + 


Hs.\-lA.n\sBA./t\s\ 


(ifa,)'  d.CIJdfcW'Pt  +  Wii)    iWifiCOil    (5-41) 


Substituting  (5.40)  into  (5.32)  produces 


d2G{q) 


=  2  Re 


Now  recall  that  asymptotically, 
and  therefore 


(£«?cr1/a  [-  (ij)*af(0Pid,(i)8j + pj*s 


-1/2 


u,  ~  7,aj 


u»  cr1/2  *  7'5f 


Of1'  u(  V         (5.42) 


(5.43) 


(5.44) 


resulting  in  ^2  being  pre-multiplied  by  ffi^P^ ,  which  is  of  course  zero.  Therefore,  (5.42) 
asymptotically  approaches 


5gg    -    2  Re  |g  |7i|2a?  [-  (a*)* df  (OPj^iJiJ  +  P**,]  a,}         (5 


=    -2  Re  J  £  |7(|2a?  (a))"  df  (1)1^3,(1)4* 


=    -2  Re 


EhHWwpyuoj 


45) 
(5.46) 
(5.47) 


since  a|a<  =  1.  To  complete  the  proof,  notice  that  the  Hessian  matrix  of  g"{q)  may  then 
be  concisely  written  as 


g"(q)  =  -2  Re 
where  D(Z)  is  defined  in  (5.29). 


£M2D"(/)Pi,D(0J 


(5.48) 
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Lemma  3   The  expected  value  E[g'(q)]  is  asymptotically  0. 
Proof: 

By  using  the  result  of  Lemma  1,  taking  the  expected  value  of  (5.14)  asymptotically 
produces 

(    L  \ 


E[g'(2)]    *    2Re{5>;D*(0Pi(,)E[w,] 
=    0 


i=i 


(5.49) 


since,  from  (4.38),  Wj  is  zero  mean.  This  may  also  be  seen  by  substituting  the  asymptotic 
value  of  u(,  7,aj,  into  (5.15)  and  observing  that  7i*a/fP±  is  zero.  I 


Lemma  4  The  expected  value  E 


g 


'(§)  (*($))' 


may  be  asymptotically  written  as 


S'(q)  (g'(?))T]  a  2  Re  {£  |7,|2LV(/)P5,D(0  j 


(5.50) 


Proof:     Consider  the  i,  j  term  of  the  matrix  g"(g).  Using  (5.27)  and  taking  the 
expected  value,  this  produces 

dGfi)dG(q) 


dqt       dtu 


E 


L      L 


(5.51) 


4  E  E  Re  {*?HmMih}  ^  {**  pi(f)a^hp} 

Using  the  identity  Re(,4)Re(B)  =  ±Re(AB  +  AB*)  [23] 

- 2  EERe E  h'-wdf wp^w^pj  a,(o 
/=i  p=i 

+7P7(dJ(P)P^(?)w;wfP^(|)dj(0]     (5.52) 
and  employing  the  conditions  on  the  interference  w  from  (4.39),  (5.52)  may  be  written  as 


E 


dG(q)dG(q) 
dq{       dqj 


E 


dG{q)dG(q) 


dq{       dqj 


Li 

2  £Re  {M2d?b)Pf,d,(0}  (5.53) 

z=i 


Equation  (5.53)  provides  the  means  to  populating  the  entire  3x3  matrix  E 
By  using  (5.29),  this  may  be  written  asymptotically  as 


E'(q)  (*(§)) 
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T~\ 


E 


g'(°q)  (g'(9))T]  a  2  Re  j  £  foW^Dtf)  j 


(5.54) 


which  concludes  the  proof. 


5.4    Bias 


The  bias  of  an  estimator  is  the  difference  between  the  expected  value  of  the  estimate 
and  the  true  value  of  the  parameter  being  estimated.  Clearly,  having  an  unknown  bias  in 
an  estimate  may  have  deleterious  impacts. 
Theorem  2   The  MLAE  attitude  estimator  is  asymptotically  unbiased,  i.e. 


E[q]-q~0 


(5.55) 


Proof:     The  bias  may  be  evaluated  by  taking  the  expected  value  of  equation  (5.11), 


since  E  [q]  —  q  =  E 


q-q 


E 


«-§) 


~    -E 


g*(Dr  g'(?) 


(5.56) 


Substituting  in  the  asymptotic  values  from  Lemma  (2)  provides 


E 


[«-?)] 


-2Re     £|7/|2D"(/)Pi(D(0 
-2Ite{EMaD"(0PiD(0 


-\  -i 


i=i 


E 


rffl) 


I'd) 


(5.57) 


(5.58) 


Now  using  the  result  from  Lemma  (3)  provides  the  following  expression  for  the  asymptotic 
bias 

-l 


[«-•> 


-2  Re     5>;|2D"(0P^D(/) 


i=i 


[0] 


(5.59) 
(5.60) 


which  shows  that  the  MLAE,  asymptotically,  is  unbiased. 


52 


5.5     Efficiency 
A  measure  of  the  quality  of  estimator  performance  is  provided  by  its  mean  squared 
error  (MSE).  The  MSE  for  the  attitude  estimate  provided  by  the  MLAE  (i.e.  the  covari- 
ance  matrix  of  the  MLAE  estimation  error),  Pq,  is  found  from 


Pq  =  E 


°\T 


(q-q)(q-q) 


(5.61) 


The  Cramer-Rao  Bound,  PCr  is  a  lower  limit  on  the  covariance  matrix  of  the  estimation 
error  for  unbiased  estimators,  and  an  estimator  that  achieves  the  Cramer-Rao  Bound  is 
said  to  be  statistically  efficient  [9].  That  is,  for  the  parameter  vector  0, 


Pi  >  PCR{e) 


(5.62) 


where  the  equality  holds  if  6  is  an  efficient  estimate  of  9. 

Theorem  3   The  MLAE  asymptotically  achieves  the  Cramer-Rao  Bound  on  estimation 

error,  and  is  therefore  statistically  efficient.   That  is, 


Pq  -  PcR{q) 


Proof:     The  estimation  error  covariance  matrix  may  be  evaluated  using  (5.11): 


(5.63) 


E 


{q-q){q-qf 


-    E 


s"(q)Y\'(q)(g'(q))T[g"(q) 


(5.64) 


since  g"(<?)  is  symmetric. 

Using  lemma  (2),  we  may  replace  the  terms  involving  g"(q)  by  their  (deterministic) 
asymptotic  values,  represented  as  g„(g): 


E 


(q-qKq-q)1 


rfS) 


&'(q)(&'Cq)j 


M) 


(5.65) 


where  (from  lemma  2) 


£{q)  =  -2Re\YJh\2i>H{l)Pip{l) 


(5.66) 
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Lemma  (4)  has  previously  addressed  the  quantity  involved  in  the  expectation,  allowing  us 
to  write 

(    L 


E[(q-q)(q-°q)T 


M) 


2  Re     £M2D"(/)Pi,D(/) 


(=i 


S'M        (5-67) 


Noting  that  the  first  two  bracketed  terms  in  (5.67)  are  inverses  (with  opposite  signs),  we 
get  the  following  expression  for  the  asymptotic  covariance  of  the  estimation  error: 


Pq    =    E 


\q-q){q-q)T] 
2Re|X>|2D"(0PiiD(0 

From  (A. 61),  the  Cramer-Rao  Bound  on  the  attitude  estimation  error  is 


(5.68) 


Pcn(q)  = 


-\  -i 


2  Re     j>|2D"(Z)Pi(;)D(0 


i=i 


(5.69) 


Observing  that  Pq  and  Pcii{q)  are  asymptotically  equal  completes  the  proof  that  the 
MLAE  is  asymptotically  efficient.  ■ 

5.6     Conclusions 

In  this  chapter  we  have  proven  that  the  Maximum  Likelihood  Attitude  Estimator 
(MLAE)  is  statically  consistent,  and  asymptotically  unbiased  and  efficient  (i.e.  achieves 
the  CRB).  These  are  three  useful  properties  that  serve  to  ensure  applicability  of  the 
algorithm  for  realistic  apphcations. 

However,  it  is  worthwhile  to  inquire  as  to  whether  an  asymptotic  (i.e.  large  sample) 
property  has  value  for  this  application.  The  simple  answer  for  GPS  attitude  determination 
is  "yes."  The  GPS  signal  strength  is  extremely  low,  so  significant  amounts  of  despreading 
and  integration  are  required  to  achieve  a  useful  SNR.  That  is,  GPS  is  simply  not  used  in  a 
small  sample  manner.  Therefore  the  asymptotic  properties  of  the  MLAE  are  useful  perfor- 
mance indicators,  since  the  expected  region  of  employment  is  a  large  sample  application. 
This  is  not  to  say  that  in  all  situations  the  MLAE  will  achieve  the  CRB,  especially  in  low 
signal  to  interference  plus  noise  (SINR)  conditions,  but  that  the  difference  between  the 
MLAE  performance  and  the  CRB  will  decrease  as  integration  time  increases. 


CHAPTER  6 
ATTITUDE  FROM  DIRECTION  FINDING 

6.1     Introduction 



In  Chapter  4  we  derived  a  jammer  resistant  maximum  likelihood  estimator  for 
the  antenna  attitude,  the  MLAE,  and  in  Chapter  5  this  estimator  was  shown  to  be 
asymptotically  unbiased,  consistent,  and  efficient.  One  drawback  to  the  MLAE  is  the 
computational  burden,  since  its  implementation  involves  a  three  dimensional  search  across 
attitude.  In  this  chapter  we  examine  methods  that,  due  to  their  computational  efficiency, 
may  prove  useful  in  a  resource  constrained  tactical  environment.  As  is  the  case  with 
suboptimal  methods,  the  computational  savings  comes  at  the  cost  of  decreased  estimator 
performance.  In  Chapter  2  attitude  determination  and  direction  finding  are  presented  as 
related  fields.  Using  this  relation,  this  chapter  develops  attitude  estimators  that  operate 
on  the  estimated  directions  to  the  satellites. 

This  chapter  is  organized  as  follows.  Section  6.2  introduces  the  general  approach  and 
form  these  estimators  take.  In  sections  6.3  and  6.4  two  formulations  of  the  estimator  are 
presented.  Finally,  section  6.5  contains  a  discussion  of  the  attributes  and  costs  of  these 
algorithms. 

6.2     Conceptual  Approach 
The  approaches  developed  in  this  chapter  employ  two  steps.  First,  the  directions 
of  arrival  of  each  satellite  are  found.  Several  methods  are  presented  in  section  2.4  that 
may  be  used  for  this  purpose.  Once  the  satellite  directions  in  the  antenna  frame  are 
estimated,  the  body  orientation  is  then  found  by  minimizing  the  error  between  the 
satellite  directions  (known  in  the  local  level  frame)  rotated  into  the  antenna  frame  and 
the  direction  estimates.  The  rotation  that  minimizes  the  error  is  chosen  as  the  antenna 
attitude. 
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To  set  up  these  direction  rinding  based  approaches,  we  define  the  following  matrices 
of  directions  where  the  subscripts  B,  and  LL  denote  the  body  frame  and  local  level  frame, 
respectively:  A  is  the  matrix  of  measured  LOS  directions  to  the  satellites,  B  is  the  known 
matrix  of  LOS  vectors  in  the  local  level  frame,  and  W  is  a  weighting  matrix. 

A    =    [£...Vl]b  (6.1) 

B    =    h-.-^JLL  (6.2) 

W    =    di&g[w1,w2,...wL]  (6.3) 

where  vt,  i  =  1, 2, . . . ,  L  is  a  3  x  1  LOS  direction  vector. 

Then  using  the  measured  and  known  directions  to  the  satellites,  the  attitude  found 
from  the  following  minimization. 

q  =  arg  min  ||  AW  -  Q(g)BW|jF  (6.4) 

This  attitude  determination  algorithm  may  be  summarized  as  follows: 
Step  1  For  each  of  the  L  satellites,  estimate  the  directions  8t,  which  correspond  to  the  LOS 

unit  vector  vLB.  Use  an  algorithm  from  Section  2.4  or  4.3. 
Step  2  Use  these  LOS  vectors  to  estimate  attitude  from  equation  (6.4). 

Two  issues  exist  in  implementing  (6.4),  once  the  directions  are  known.  First,  the 
equation  is  non-linear  in  the  attitude  q,  and  therefore  the  attitude  cannot  be  found 
directly  through  matrix  inversion  or  psuedo-inversion.  Second,  the  choice  for  the  weighting 
matrix  is  undefined.  We  address  the  first  issue  through  an  iterative  method  based  on 
a  linearized  version  of  (6.4).  To  address  the  second  issue,  two  options  of  the  weighting 
matrix  are  presented.  For  the  remainder  of  this  chapter  we  parameterize  attitude  only  as 
a  unit-norm  quaternion  to  avoid  the  convergence  issues  and  singularities  that  are  possible 
when  using  Euler  angles. 
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6.3     Equal  Satellite  Weighting 

As  in  Markel  et  al.  [4],  one  choice  for  the  weighting  matrix  is  the  identity  matrix. 

With  this  choice  every  satellite  contributes  equally  to  the  solution,  and  (6.4)  may  be 

written  as 

g  =  argmin||A-Q(?)B||F  (6.5) 

i 

Conceptually,  the  true  LOS  vectors  in  body  frame  may  be  considered  to  have  been  rotated 
by  the  true  attitude  q.  Using  this  concept  we  define  Z(q)  as  the  matrix  of  LOS  vectors  in 
the  local  level  frame  rotated  by  the  true  attitude,  and  Z(q)  as  the  matrix  of  LOS  vectors 
in  the  local  level  frame  rotated  by  the  attitude  represented  by  the  arbitrary  quaternion  q. 

Z(q)  =  Q(q)B  (6.6) 

Z(q)  =  Q(q)B  (6.7) 

Now  A,  the  matrix  of  measured  LOS  vectors  in  body  frame,  may  be  considered  to 
be  an  estimate  of  Z(q),  say  Z(q),  since  the  directions  are  measured  in  the  body  frame 
and  known  in  the  local  level  frame.  In  order  to  extract  the  attitude  from  the  measured 
directions  in  A,  i.e.  to  estimate  q  from  the  estimate  of  Z(q),  substitute  Z(q)  (i.e.  A)  and 
Z(q)  into  (6.5). 


q  =  arg  mm 


[Z(q)  -  Z(q) 


(6-8) 


In  the  no-noise  case  it  is  easy  to  see  that  the  quaternion  that  minimizes  (6.8)  is  q  =  q, 
since  Z(q)  will  equal  Z(q). 

This  quaternion  q  of  (6.8)  may  be  found  through  an  iterative  process.  This  procedure 
follows  the  general  form  of  Cohen  [7]  (which  was  solving  for  attitude  in  Euler  angles 
by  using  phase  differences),  but  has  been  rewritten  for  this  application.  We  begin  by 
expressing  Z(q)  through  its  Taylor  series  expansion. 

Z(q)  =  Z(q)  +  J2  -S^*  (q  -  q)  +  {higher  order  terms}  (6.9) 


Ignoring  the  higher  order  terms,  (6.9)  may  be  written  as 


where 


E  =  X8q 


E    =    vec(Z(q)  -  Z(q) 
X    = 


.dZ(q).  ,0Z(q),  ,SZ(q),  ,0Z(q) 

vec(    a_     ),  vec(    Q       ),  vec(    a       ),  vec( 


<9<7i 


<9<72 


dq3 


dq4 


6q  q-q 

Sq, 

£q2 
£q3 
£q4 
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(6.10) 


(6.11) 
(6.12) 
(6.13) 

(6.14) 


and  vec(-)  represents  stacking  all  columns  of  a  matrix  into  a  single  column  vector.  The 
quaternion  update  vector  is  found  from  the  psuedo-inverse  of  X. 


<Jq=(X*X)   JX"E  (6.15) 

The  iterative  algorithm  for  solving  (6.8)  may  be  summarized  as  follows: 
Step  1  Choose  an  initial  value  for  the  estimated  attitude  quaternion  q.  This  may  be  the 
last  attitude  estimated  from  the  previous  data  set,  or  if  no  a  priori  information 
exists,  the  unit  quaternion. 
Step  2  Calculate  the  quaternion  update  £q  from  (6.15). 
Step  3  Calculate  the  new  quaternion  by  adding  q  to  6q. 
Step  4  Normalize  the  new  quaternion  to  unit  norm. 

Step  5  Return  to  step  2  until  convergence  occurs,  i.e.  when  8q  is  arbitrarily  small. 
It  is  important  to  include  Step  4,  since  the  above  method  does  not  intrinsically  ensure 
that  the  quaternion  remain  unit  norm.  This  algorithm  will  be  referred  to  in  this  work  as 
DFU,  since  it  uses  direction  finding  to  obtain  the  attitude  estimate,  and  does  not  use  an 
weighting  matrix  that  is  the  identity  matrix  (i.e.  is  unweighted). 
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6.4    Satellite  Weighting  via  Adapted  SINR 
The  algorithm  developed  in  this  section  uses  the  same  iterative  approach  as  the 
DFU  algorithm,  however  now  the  weighting  matrix  is  not  the  identity  matrix.  The  choice 
employed  here  is  to  weight  each  satellite's  contribution  to  the  attitude  estimate  by  its 
adapted  SINR.  The  adapted  SINR  is  found  from  the  array  response  vector  to  the  source, 
the  interference  covariance  matrix,  and  the  complex  gain  [42] . 1 

SmRAdapted  =  171^(0)0-^(0)  (6.16) 

where  7  is  the  complex  gain,  a(0)  is  the  array  response  vector  from  the  source  located  at 
direction  9,  and  C  is  the  interference  covariance  matrix.  In  practice,  these  values  may  be 
replaced  with  their  estimates. 

Throughout  this  chapter  the  method  of  determining  the  directions  to  the  satellites 
in  the  antenna  frame  has  been  left  unspecified.  However,  if  the  maximum  likelihood 
algorithms  of  (4.28)  and  (4.29)  of  Section  4.3  are  employed  to  estimate  the  directions 
to  the  satellites,  then  calculation  of  the  adaptive  SINR  is  essentially  provided  by  the 
algorithm  already.  Notice  that  the  adaptive  SINR  can  be  viewed  in  terms  of  the  whitened 
array  response  vector  a. 

a  =  cHa  (6.17) 

and  therefore 

SWRAdapted  =  |7|2|a|2  (6.18) 

The  greater  the  attenuation  from  the  whitening,  the  lower  the  adapted  SINR.  We  will 
refer  to  this  algorithm  as  DFW,  since  it  uses  direction  finding  and  a  weighting  matrix. 


1  The  term  adapted  SINR  is  actually  not  used  in  Applebaum  [42],  but  is  a  common 
term  today  used  to  represent  the  signal  to  interference  plus  noise  ratio  after  application  of 
the  optimal  beamformer  weights  derived  in  Applebaum  [42] . 
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6.5    Conclusions 
Three  facets  of  this  method  of  attitude  estimation  warrant  further  discussion.  The 
first  is  the  relationship  among  the  satellites.  In  a  multi-user  wireless  communications 
system  where  the  individual  sources  are  essentially  unconstrained  in  their  possible 
location,  determining  the  directions  to  each  source  separately  is  acceptable,  and  in 
most  cases  desirable  as  it  reduces  the  dimensionality  of  the  search.  However,  the  GPS 
source  (satellite)  locations  are  constrained  in  location  by  their  known  flight  paths.  The 
constellation,  in  the  local-level  frame,  is  known  by  the  user  receiver,  as  is  the  receiver 
position,  providing  additional  information  that  could  be  used  in  the  direction  estimation 
process.  Indeed,  the  direction  to  one  satellite  is  not  completely  independent  from  the 
directions  to  all  others.  Therefore,  to  estimate  the  directions  independently,  as  in  (4.29),  is 
a  computationally  appealing  but  suboptimal  method  of  attitude  estimation.  In  Chapter  8, 
simulation  results  quantify  the  degradation  of  performance  of  these  suboptimal  algorithms 
from  the  optimal  MLAE. 

The  second  area  of  discussion  centers  around  the  contribution  of  each  satellite.  In 
section  4.6  a  graphical  example  was  used  to  show  that  the  contribution  of  each  satellite 
is  not  the  same  across  attitude-space.  Inherent  in  the  metric  of  equation  (4.59)  is  the 
(optimal)  mechanization  of  how  each  satellite  contributes  across  attitude  space  to  the 
attitude  estimate,  even  in  the  presence  of  interference  and  jammers.  With  the  LOS 
mapping  approach,  this  optimality  is  lost.  Using  equal  weighting,  each  satellite  not  only 
contributes  equally  across  attitude-space,  they  all  contribute  equally  relative  to  each  other. 
With  the  adapted  SINR  weighting  satellites  near  the  jammer  sources  will  not  be  allowed 
to  contribute  as  much  to  the  attitude  estimate  as  those  with  larger  angular  separations. 
However,  the  contribution  from  any  satellite  is  still  equal  across  attitude-space,  since  only 
a  scalar  weighting  is  apphed  to  the  satellite.  One  can  clearly  see  how  this  is  a  suboptimal 
approach  by  trying  to  represent  the  information  of  Figures  4.3,  4.4,  and  4.5  with  a  single 
scalar  per  figure. 

The  third  discussion  topic  is  the  relative  performance  of  these  two  new  attitude 
estimators.  In  general,  weighting  the  satellite  contributions  by  the  adapted  SINR  in 
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the  attitude  minimization  process  outperforms  equal  satellite  weighting.  As  shown  in 
Chapter  8,  this  is  often  the  case  for  the  scenarios  simulated.  However,  in  some  of  the 
cases  examined  the  weighted  estimator  has  a  slightly  larger  mean  total  error  than  the 
unweighted.  This  seems  to  occur  when  the  adapted  SINR's  are  below  approximately  6-10 
dB,  which  occurs  at  the  fastest  update  rates  (i.e.  smallest  integration  times),  with  the 
fewest  sensors,  and  smallest  baseline  lengths.  It  seems  that  in  this  case,  the  weighting  puts 
too  great  a  confidence  in  the  noisy  measurements  of  a  few  satellite  directions. 

Although  not  investigated  in  this  work,  one  possible  improvement  to  this  process 
would  be  to  use  a  switchable  weighting  scheme  based  on  the  estimated  adapted  SINR.  As 
discussed  in  Chapter  9,  this  scheme  would  use  the  weights  only  when  the  adapted  SINR's 
crossed  an  empiricaUy  derived  threshold. 


CHAPTER  7 
WIDER  BASELINES  AND  DUAL  FREQUENCY  USE 

7.1     Introduction 

It  is  well  known  that  as  the  separation  between  sensors  increases  the  variance  of  the 
direction  estimation  error  decreases?    However,  this  increase  in  sensor  spacing  comes  at 
a  cost.  If  the  sensors  are  spaced  farther  than  the  Nyquist  criterion,  grating  lobes  appear, 
resulting  in  an  ambiguity  of  direction. 

A  similar  phenomena  arises  for  the  attitude  determination  application.  With  the 
conventional  AD  approach  of  using  phase  differences,  only  the  fractional  portion  of  the 
total  phase  may  be  measured.  However,  the  total  phase  difference  between  any  two 
sensors  receiving  a  satellite  signal  is  composed  of  the  fractional  part  (i.e.  the  measured 
portion),  and  the  integer  portion.  Since  the  integer  portion  is  not  measured,  yet  must  be 
accounted  for  to  properly  determine  the  antenna  attitude,  it  is  referred  to  as  the  integer 
ambiguity  by  the  attitude  determination  community. 

Depending  on  the  sensor  spacing  and  the  number  of  satellites  used  to  determine 
attitude,  the  integer  ambiguity  may  be  resolved  uniquely.  However,  this  typically  requires 
an  intensive  search  of  all  the  possible  ambiguities.  One  method  to  eliminate  many  of  the 
ambiguities  is  through  a  process  called  "widelaning."  The  widelane  phase  measurement  is 
formed  from  the  difference  of  the  phase  measurements  from  the  two  GPS  frequencies,  i.e. 
the  LI  and  L2  frequencies: 

^widelane  =  $L1  -  $L2  (7.1) 


1  This  fact  is  easily  demonstrated  by  the  CRB  of  the  direction  estimate  vx  for  a  two 
sensor  array:  o2Ux  >  1/[SNR  n2o9},  where  d  is  the  distance  between  the  sensors. 
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This  gives  an  effective  wavelength  of  [37] 

^widelane  —  T  T  (7  .z) 

which  for  the  two  GPS  frequencies  of  1.228  and  1.575  GHz  corresponds  to  .86  meters,  or 
approximately  four  times  the  wavelength  corresponding  to  the  highest  GPS  frequency. 
This  significantly  reduces  the  size  of  the  integer  ambiguity  search.  The  price  paid  for  the 
increased  effective  wavelength  is  a  doubling  of  the  receiver  noise,  since  the  receiver  noise  in 
each  frequency  channel  is  uncorrelated. 

The  Maximum  Likelihood  Attitude  Estimator  (MLAE)  of  Chapter  4  is,  in  many 
cases,  unaffected  by  the  integer  ambiguity  phenomenon.  Since  the  MLAE  searches 
attitude  with  a  fixed  array  configuration,  it  can  often  perform  unambiguously  even  when 
the  sensors  are  farther  apart  than  the  Nyquist  criterion.  Provided  that  the  uncertainty 
in  attitude  over  which  to  search  is  not  too  large,  the  MLAE  will  converge  to  a  single 
local  minima  within  the  uncertainty  when  several  satellites  are  being  tracked,  even  if  the 
sensor  spacing  is  greater  than  the  Nyquist  criterion.  This  occurs  because  inclusion  of  the 
additional  satellites  (over  the  minimum  two  required)  causes  the  MLAE  metric  at  the 
bogus  attitudes  to  indicate  that  these  attitudes  are  less  likely  than  the  true  attitude.  If 
the  MLAE  is  cast  as  a  minimization  as  in  equation  (4.49),  the  metric  value  at  the  bogus 
attitudes  will  not  typically  be  as  small  as  the  value  at  the  true  attitude. 

However,  a  situation  similar  to  the  integer  ambiguity  problem  does  exist  for  the 
MLAE  when  the  number  of  satellites  being  tracked  is  small,  the  sensor  spacing  is  greater 
than  the  Nyquist  criterion,  and  the  attitude  uncertainty  (i.e.  the  region  to  search)  is 
large.  In  this  case  the  MLAE  metric  at  several  regions  of  attitude  space  may  be  close 
to  the  value  at  the  true  attitude,  resulting  in  local  minima  of  (4.49)  that  could  result  in 
false  solutions  to  a  minimization  search.  Therefore,  it  would  be  desirable  to  develop  a 
dual  frequency  version  of  the  MLAE  developed  in  chapter  4  that  reduces  the  number  of 
false  attitude  solutions  like  the  conventional  attitude  determination  process  of  widelaning. 
Unfortunately,  a  simple  "phase  difference"  approach  like  widelaning  is  not  applicable  to 
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the  MLAE.  Since  the  MLAE  does  not  operate  on  phases  alone  but  on  the  complex  values 
of  the  data  vectors,  a  more  elaborate  estimation  process  is  required. 

In  this  chapter  we  will  develop  a  method  of  incorporating  the  second  GPS  frequency 
into  the  MLAE.  We  will  show  that  this  dual  frequency  extension  provides  two  positive 
attributes:  a  decrease  in  estimation  error  and  a  reduction  in  possible  false  attitudes. 
The  following  section  derives  the  dual  frequency  MLAE.  Similar  to  the  initial  MLAE 
derivation,  this  begins  with  a  new  look  at  the  array  response  vector.  Following  the 
derivation  we  present  a  section  that  examines  the  reduction  in  false  attitudes  by  using 
the  second  GPS  frequency.  Finally,  we  show  that  this  dual  frequency  attitude  estimator 
asymptotically  achieves  the  Cramer-Rao  Bound,  and  produces  higher  quality  estimates 
than  the  single  frequency  version. 

7.2     Dual  Frequency  Maximum  Likelihood  Attitude  Determination 
Before  this  section  develops  the  dual  frequency  attitude  estimator,  first  review  the 
single  frequency  attitude  estimator  of  Chapter  4.  This  estimator  used  the  L  demodulated 
data  vectors  (each  M  x  1),  the  L  interference  covariance  matrices,  and  an  array  response 
vector  parameterized  on  the  antenna  attitude  (and  using  the  known  directions  to  the 
satellites  in  the  local  level  frame)  to  define  the  log  likelihood  surface.  The  minima  or 
maxima  (depending  on  whether  equation  (4.49)  or  (4.51)  are  used)  of  this  function  is 
chosen  to  be  the  attitude  estimate. 

Although  the  development  of  this  estimator  is  performed  using  baseband  signals, 
the  carrier  frequency  of  the  GPS  signals  appears  in  the  array  response  vector  of  equation 
(4.31)  through  the  division  by  the  wavelength  A.  Assuming  that  the  conditions  on  the 
interference  from  Chapter  4  are  met  for  both  GPS  frequencies,  the  derivation  of  this 
estimator  will  apply  equally  well  to  either  frequency,  since  the  statistical  properties  of  the 
DS-SS  satellite  waveforms  are  the  same  for  both  the  LI  and  L2  GPS  frequency.  Without 
introducing  any  additional  unknowns  into  the  estimator,  the  array  response  vector  may  be 
written  as  a  function  of  the  unknown  attitude  q,  the  known  satellite  direction  in  the  local 
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level  frame  vll,  and  now  the  known  wavelength  A: 


a(^LL,  A,  q)  = 


(7.3) 


Now  consider  a  situation  that  makes  use  of  this  modified  array  response  vector.  A 
receiver  capable  of  demodulating  both  frequencies  of  the  GPS  P(Y)  code  can  generate  2L 
data  vectors,  L  per  frequency  for  each  integration  period.  To  distinguish  the  data  vectors 
gathered  from  the  different  frequencies,  the  notation  representing  the  data  vector  obtained 
from  demodulating  with  the  spreading  sequence  from  satellite  /  is  expanded  from  u(  to 
u(',/)'  where  /  denotes  which  of  the  two  possible  carrier  frequencies.  Using  this  expanded 
notation,  we  define  the  augmented  space-satellite  data  vector,  Ua,  of  the  received  data 
similar  to  equation  (4.40): 


UQ  = 


(7.4) 


u(i,i) 

U(2,l) 

U(£,l) 
"(1,2) 
U(2,2) 

"(L,2) 

Letting  U,  represent  the  space-satellite  data  vector  for  frequency  i,  i  =  1, 2,  equation  (7.4) 
may  be  written  in  block  form  as 

u, 


ua  = 


u, 


(7.5) 


Now  consider  the  statistics  of  the  interference.  Similar  to  equation  (7.5),  the  aug- 
mented interference  vector  (containing  thermal  noise,  jammers,  and  the  multiple  access 
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interference  from  other  GPS  satellites,  as  in  Chapter  4)  may  be  written  in  block  form  as 


W„ 


Wi 
W2 


(7.6) 


where  W<  represents  the  ML  x  1  interference  vector  for  frequency  i,  i  =  1,2.  Assuming 
that  the  jammer  waveforms  are  uncorrelated  frequency  to  frequency,  Rss_a,  the  covari- 
ance  of  W„,  is  block  diagonal  just  as  KS3,  the  covariance  of  W  for  a  single  frequency 
in  equation  (4.46)  is  block  diagonal.  This  implies  that  the  entire  interference  at  one 
frequency/satellite  channel  is  uncorrelated  from  the  interference  in  all  other  frequen- 
cies/satellite channels.  Letting  Cij  represent  the  covariance  of  the  interference  in  the  Ith 
satellite  channel,  I  —  1,2, . . .  ,L  and  the  /th  frequency,  i  =  1,2,  the  augmented  interference 
covariance  matrix  may  be  written  as 


R, 


Rss_i        0 
0       Rss-2 

cu     0  0 

0      C2,i      0  0 

o      •••    ••.  0 

o        •••      0  d,i 


(7.7) 


Cll2         0         •••  0 

0      Co,     0        0 


'2,2 


0 
0 


(7.8) 


••      '•.       0 
...        0      CL,2 

Using  the  above  definitions,  the  asymptotic  log  likelihood  ratio  of  (4.47)  may  be  written  in 
dual-frequency  form  as 


i=i 


^R-EKy)"^)^>A>'«)]*cS)[u('.i)-Tf(y)«('i)Ai,?)]       + 

L 

]£Hw  -  7(/,2)a(i/(,  A2,  q)]HC-^[u{l}2)  -  7,i2a(n,  A2,  q)]    (7.9) 
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As  with  the  single  frequency  attitude  estimator,  a  closed  form  expression  for  7^/)  at  a 
stationary  point  may  be  obtained  by  setting  the  partial  derivative  of  (7.9)  with  respect  to 
7((j)  equal  to  zero  and  solving  for  7^/)  in  terms  of  q. 

a//(^,A/,g)C^1/)u(iJ) 


7(i,/) 


(7.10) 


Substituting  in  the  complex  gain  expression  of  (7.10)  provides  the  following  minimiza- 
tion for  the  attitude  q: 


/=i 


q  =  argmin^[u(U)  -  7(U)*K  Kq^C^u^  -  %tl)a(v,,  Xlfq)]         + 

L 

X][u(',2)  -  7((,2)a(^,  A2,g)]//C(712)[u(,i2)  -  7,,2a(^,  X2,q)}    (7.11) 


/=i 


By  expanding  and  removing  terms  that  do  not  vary  with  attitude,  and  substituting 
in  the  expression  for  7(ii/)  from  equation  (7.10),  equation  (7.11)  may  be  written  as  a 
maximization: 


q  =  arg  max  } 

a       ^— ' 


rtt  I  la"KAi><?)c(u)aKAi,g) 


+ 


\&{vi^2,q)}"C(l]2)u(la)\2 
af/(^,A2,  q)C^)a(ui,X2,q) 


(7.12) 


As  with  the  single  frequency  MLAE,  a  simpler  form  of  the  expression  may  be  written 
using  the  shortened  form  of  the  array  response  vector: 


a(',/)(g)  =  aKA/,g) 


(7.13) 


And  as  in  (4.53)  we  may  define  the  "whitened"  dual-frequency  attitude  array  response 
vector: 


hi  j)  =  c((,1/)2aK'V,<7) 


(7.14) 
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where 


'('./) 


_      p.l/2  p.l/2 


pl/2 


Using  these  definitions,  (7.12)  may  be  rewritten  as 

l</)CJ)2a('./)(<?)|: 


(7.15) 
(7.16) 


L       2 

arg  max  Y^  V^ 

9     i=i  /=i 
z-      2 


N,/)(?)l2 

""?*  £  S  U^/)C(i)1//)2p«(',/)(9)C(/,/)2u('./) 


/=1   /=1 


where 


ac,/) 


(9) 


_   *(',/)  (<?)  5"/)(g) 

5?,/)(9)  »(<,/)  (?) 


(7.17) 
(7.18) 

(7.19) 


The  attitude  estimators  of  equations  (7.11),  (7.12),  and  (7.31)  will  be  referred  to  in  this 
work  as  the  dual- frequency  MLAE. 

Several  observations  may  be  made  about  the  dual-frequency  MLAE.  First,  as  with  the 
single  frequency  version,  an  expression  for  the  complex  gain  may  be  found  using  data  from 
a  single  satellite  and  single  frequency.  However,  the  attitude  estimate  includes  information 
from  all  satellites  and  both  frequencies.  Although  coupled,  this  could  be  implemented 
in  a  parallel  structure,  having  perhaps  one  computer  calculating  the  metric  for  a  set  of 
attitudes  using  the  data  from  the  first  frequency  while  another  processes  the  data  from 
the  second  frequency.  The  metric  values  could  then  be  combined  for  the  set  of  attitudes 
investigated,  and  the  attitude  with  the  minimum  combined  metric  value  would  be  chosen 
as  the  estimate. 

Although  the  metric  was  derived  using  the  same  satellites  on  both  bands,  this  in  not 
a  requirement.  In  fact,  none  of  the  satellites  need  to  appear  at  both  frequencies  for  the 
algorithm  to  operate.  Similarly,  this  approach  performs  if  the  jammer  scenario  is  different 
on  the  two  frequencies.  For  example,  if  one  band  was  jammed  and  the  other  unjammed, 
this  algorithm  would  perform  better  than  a  single  frequency  MLAE  operating  at  either  the 
jammed  or  even  the  unjammed  frequency. 
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7.3     Reduction  in  False  Attitudes 
Consider  the  situation  where  only  a  few  satellites  are  visible,  a  large  uncertainty 
in  attitude  exists,  and  to  achieve  higher  resolution  attitude  estimates  the  sensors  are 
spaced  farther  apart  than  a  half  wavelength.  In  practice,  this  scenario  may  occur  either 
on  startup  or  if  INS  aiding  is  either  not  available  or  not  operating.  In  this  case,  there  may 
be  several  areas  of  attitude  space  that  lead  to  local  minima  of  equation  (7.11)  (or  maxima 
of  equations  (7.12)  and  (7.31)).  To  find  the  attitude  estimate,  the  extrema  of  each  area 
would  have  to  be  found  and  sorted,  with  the  attitude  estimate  chosen  as  the  attitude 
corresponding  to  the  global  extrema.  It  is  clear  that  a  reduction  in  the  number  of  areas 
that  may  contain  an  extrema  (perhaps  crossing  a  threshold  and  therefore  requiring  that 
the  extrema  be  found  and  its  metric  value  obtained  for  comparison  to  others)  would  be 
desirable. 

To  visualize  how  incorporation  of  the  second  frequency  reduces  the  number  of  false 
attitude  positions,  consider  the  simpler  direction  finding  application.  Figure  7.1  shows  the 
(normalized)  likelihood  ratio  for  a  received  signal  impinging  on  the  array  from  broadside 
for  a  15  element  uniform  linear  array  (ULA)  with  a  2A  spacing  between  elements.  In  this 
example  no  jammers  exist,  i.e.  the  likelihood  value  is  proportional  to  the  quiescent  array 
pattern.  Notice  that  several  directions  produce  likelihood  values  equal  to  the  maximum  at 
the  correct  direction.  For  a  radar  these  grating  lobes  cause  two  deleterious  effects:  they 
produce  an  ambiguity  in  the  direction  estimation  process,  and  they  allow  undesired  signals 
from  these  ambiguous  directions  to  pass  through  the  antenna  with  the  same  gain  as  the 
desired  signal.  Now  consider  the  likelhood  value  for  the  same  array  receiving  signals  at 
a  carrier  frequency  equal  to  1.25  times  greater.  From  Figure  7.2  it  can  be  seen  that  the 
locations  of  the  grating  lobes  have  changed.  The  extension  of  this  example  to  the  dual- 
frequency  MLAE  is  obtained  by  summing  the  likelihood  values  obtained  from  these  two 
frequencies,  as  shown  in  Figure  7.3  .  In  this  case  the  grating  lobes  are  attenuated,  since 
only  the  true  direction  appears  at  the  same  angular  location  as  frequency  changes. 


Array  Response  for  Frequency  1 
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Figure  7.1:  Normalized  likelihood  value  vs.  direction  of  arrival  for  a  15  element  ULA  re- 
ceiving a  narrowband  signal  at  a  frequency  corresponding  to  a  2A  spacing.  True  direction 
is  on  boresight,  i.e.  0  degrees. 


Array  Response  for  Frequency  2 
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Figure  7.2:  Normalized  likelihood  value  vs.  direction  of  arrival  for  a  15  element  ULA 
receiving  a  narrowband  signal  at  a  frequency  corresponding  to  a  ^A  spacing.  True  direc- 
tion is  on  boresight,  i.e.  0  degrees. 
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Combined  Array  Response  for  Both  Frequencies 
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Figure  7.3:  Likelihood  value  vs.  direction  of  arrival  for  a  15  element  ULA  incorporating 
data  collected  at  two  frequencies  corresponding  to  2A  and  j^A  spacing.  True  direction 
is  on  boresight,  i.e.  0  degrees.  The  grating  lobes  are  reduced  in  amplitude  from  the  true 
direction  of  arrival,  producing  a  clear  indication  which  direction  is  correct. 
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Now  consider  how  this  applies  to  the  attitude  determination  case.  In  the  following 
example,  the  true  attitude  (in  Euler  angles)  is  [0,0,0],  and  the  attitude  uncertainty  is 
assumed  to  be  ±88  degrees  in  each  coordinate.  Three  likelihood  surfaces  are  generated 
at  a  two  degree  resolution:  the  MLAE  (equation  (4.49))  using  frequency  LI  ,  the  MLAE 
(equation  (4.49))  using  frequency  L2,  and  the  dual  frequency  MLAE  (7.11).  Each  of 
these  estimators  chooses  the  attitude  that  minimizes  the  likelihood  expression,  and  each 
likelihood  surface  is  calculated  using  two  satellites.  A  threshold  of  lOdB  below  the  mean 
value  of  each  surface  was  chosen,  and  all  points  in  attitude  space  below  this  threshold  were 
considered  possible  candidates  for  the  attitude  estimate.  In  order  to  show  the  reduction  in 
false  values  of  attitude,  any  point  that  had  a  contiguous  connection  to  the  true  attitude 
was  removed  in  the  following  plots.  Any  point  with  such  a  connection  would  lead,  by  a 
simple  minimization  scheme,  to  the  true  attitude.  Therefore,  for  comparison  purposes  we 
wish  to  examine  only  those  points  that  would  lead  to  false  values  of  attitude,  i.e.  a  local 
minimum  that  is  not  at  the  true  attitude.  The  intent  of  this  is  to  evaluate  the  potential 
decrease  in  false  attitudes  when  using  the  dual-frequency  MLAE  over  using  a  single 
frequency  MLAE. 

Figure  7.4  shows  a  three  dimensional  representation  (over  the  three  Euler  coordinates 
of  attitude)  of  the  attitudes  where  the  likelihood  value  was  below  the  threshold  using 
the  LI  frequency.  To  provide  an  easier  visualization  of  these  points,  Figure  7.5  shows 
a  two  dimensional  version  of  the  same  plot,  obtained  by  projecting  all  points  onto  the 
it,  =  0  plane.  As  discussed  above,  the  attitude  values  that  connect  to  the  true  attitude 
are  removed,  leaving  only  points  that  would  lead  to  a  a  false  attitude.  To  illustrate  this, 
figure  7.5  shows  all  attitudes  with  metric  values  below  the  threshold,  and  figure  7.6  shows 
only  those  that  do  not  connect  to  the  true  attitude.  The  relatively  few  points  leading  to 
the  true  attitude  have  been  "thinned  out."  The  jammer  scenario  used  is  the  same  as  that 
involving  the  first  jammer  in  the  random  jammer  study  of  Chapter  8.  The  four  sensors 
used  are  spaced  in  a  "quad"  formation  with  spacing  of  approximately  .6  meters. 

Figure  7.7  presents  the  false  attitudes  for  the  MLAE  using  the  L2  frequency  in  the 
same  manner.  In  Figure  7.8,  which  contains  the  likelihood  points  below  the  threshold 
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for  the  dual-frequency  MLAE,  it  is  apparent  that  a  significant  reduction  in  the  number 
of  points  leading  to  false  attitude  solutions  exists.  Similar  to  the  grating  lobe  reduction 
in  the  direction  finding  application  above,  the  dual-frequency  MLAE  "smoothes"  the 
likelihood  value  away  from  the  true  value  of  attitude,  reducing  the  number  of  points  with 
likelihood  values  that  may  be  considered  candidates  for  the  global  minimum. 

To  illustrate  the  reduction  in  false  attitudes,  three  additional  data  sets  are  presented. 
Figures  7.9,  7.10,  and  7.11  present  the  single  frequency  MLAE  for  LI  and  L2,  and  the 
dual-frequency  MLAE  using  satellites  1  and  5.  Similarly,  figures  7.12,  7.13,  and  7.14 
present  data  using  satellites  1  and  6,  and  figures  7.15,  7.16,  and  7.17  present  data  using 
satellites  2  and  3.  Table  7.1  contains  the  number  of  false  attitude  possibilities  for  the  three 
estimators  for  several  more  satellite  combinations. 

Raw  Points  for  Satellites  1  and  4.  Frequency  1 
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Figure  7.4:  Three  dimensional  plot  of  all  attitudes  whose  metric  value  is  below  the  lOdB 
threshold,  using  Satellites  1  and  4,  and  frequency  LI. 
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Raw  Points  for  Satellites  1  and  4.  Frequency  1 
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Figure  7.5:  Two  dimensional  view  of  all  attitudes  whose  metric  value  is  below  the  lOdB 
threshold,  using  Satellites  1  and  4,  and  frequency  LI. 
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Figure  7.6:  Two  dimensional  view  of  those  attitudes  whose  metric  value  is  below  the  lOdB 
threshold  that  are  not  contiguous  to  the  true  attitude,  (i.e.  false  solutions)  using  Satellites 
1  and  4,  and  frequency  Ll. 
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Satellites  1  and  4.  Frequency  2 
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Figure  7.7:  False  attitude  solutions  using  Satellites  1  and  4,  and  frequency  L2. 
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Figure  7.8:  False  attitude  solutions  using  Satellites  1  and  4,  and  frequency  the  dual  fre- 
quency MLAE 
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Satellites  1  and  5.  Frequency  1 

mpr       nil    ||«{ji 

IP1     ***   ir     jL, 

:;iii::.     I  Hi'*/   !«Vl* 


u 


l<       ::j: 

'    jl        :: 


ti 


Ha  Mj 


::: 
iih- 


1  Hi  v1'  i 

tf.JLfc...] 


1     i 


-80        -60        -40        -20  0  20  40  60  80 

Roll  (Degrees) 

Figure  7.9:  False  attitude  solutions  using  Satellites  1  and  5,  and  frequency  LI. 
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Figure  7.10:  False  attitude  solutions  using  Satellites  1  and  5,  and  frequency  L2. 
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Figure  7.11:  False  attitude  solutions  using  Satellites  1  and  5,  and  the  dual  frequency 
MLAE. 
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Figure  7.12:  False  attitude  solutions  using  Satellites  1  and  6,  and  frequency  LI. 
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Figure  7.13:  False  attitude  solutions  using  Satellites  1  and  5,  and  frequency  L2. 
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Figure  7.14:  False  attitude  solutions  using  Satellites  1  and  5,  and  the  dual  frequency 
MLAE. 
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Figure  7.15:  False  attitude  solutions  using  Satellites  2  and  3,  and  frequency  LI. 
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Figure  7.16:  False  attitude  solutions  using  Satellites  2  and  3,  and  frequency  L2. 
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Satellites  2  and  3.  Both  Frequencies 
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Figure  7.17:  False  attitude  solutions  using  Satellites  2  and  3,  and  the  dual  frequency 

MLAE. 


Table  7.1:  Number  of  false  solutions,  i.e.  points  below  the  "possible  minimum"  threshold 
that  do  not  converge  to  the  true  solution.  Resolution  in  each  Euler  angle  is  two  degrees. 
Scenario  is  the  same  as  that  involving  the  first  jammer  in  the  random  jammer  study  of 
Chapter  8. 


Satellite  Combination 

MLAE  LI 

MLAE  L2 

Dual  Freq.  MLAE 

1-4 

3681 

5052 

683 

1-5 
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4690 

153 
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7.4    Dual-Frequency  MLAE  Performance  Increase 
In  this  section  we  analytically  examine  the  performance  of  the  dual-frequency  MLAE. 
Specifically,  we  will  show  that  the  dual-frequency  estimator  may  be  written  in  the  form 
similar  to  that  used  for  the  single  frequency  estimator.  Then,  using  the  results  of  Chapter 
5,  we  show  that  the  dual-frequency  MLAE  has  the  same  desirable  properties  as  the  single 
frequency  MLAE. 

We  begin  be  examining  the  dual-frequency  MLAE  from  equation  (7.11),  repeated  here 
for  convenience. 


q  =  argmin  £[u(U)  -  7(U)a(i/;,  Ai.^C^foty,,  -  7(U)*K  Ai,  ?)]         + 


^[u((,2)  -  7«,2)a(i/,,  A2, 9)]i/C((12)[u(i,2)  -  7,,2a(^,  A2,  q)}    (7.20) 
This  equation  may  be  rewritten  in  simpler  notation  as 

2L 

q  =  argmin  J}^  -  fy&cMf&ffa  -  7cM<7)]  (7.21) 

where  we  have  defined 


"    C-l 
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u<  =  < 
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a("<>Ai>9)  1<C<-^ 


a(i/c_L,A2t?)    L  +  1<C<2L 

"(CD  1  <  C  <  L 

u(C-£,2)    L  +  1  <  C  <  2L 


(7.22) 


(7.23) 


7(c,D  1  <  C  <  L 

TC  =  \  (7.24) 

7(c-l,2)    I  +  1  <  C  <  2L 


C(Cfi)         1  <  C  <  i 
C((-£,i)    I  +  1  <  C  <  2L 
With  this  substitution,  the  frequency  is  no  longer  an  argument  of  the  array  response 
vector,  the  data  vector,  or  the  covariance  matrices.  The  terms  relating  to  the  second 


(7.25) 
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frequency  are  simply  indexed  as  the  last  L  C's.  With  this  notation,  we  can  rewrite  the 
maximization  equation  (7.12)  as 


q  =  arg  max 


2L 

£ 

l=i  L 


\*dq)HC7%\* 


^(q)BC7%(q) 


(7.26) 


As  in  Chapter  4,  we  can  simplify  further  by  defining  the  whitened"  dual-frequency 
attitude  array  response  vector: 

4(f)  =  ^%(q)  (7.27) 

where 


c<   =  6-"26-»2 

Using  these  definitions,  (7.12)  may  be  rewritten  as 


C"1/2 


(7.28) 
(7.29) 


2L    ixHg-l/a*. 


EIUC  ^c 
I4(?)l2 


4fo)lS 


=    max 


C-l 

2L 


1/2 


X-l/2v 


where 


<M<7)  —   X// 


4(g)  gw 

951  4(9) 


(7.30) 
(7.31) 

(7.32) 


Notice  that  equation  (7.31)  is  in  exactly  the  same  form  as  the  single  frequency  version 
of  equation  (4.56).  With  this  form  we  can  use  the  same  asymptotic  formulas  for  efficiency 
as  in  Chapter  5  and  Appendix  A.  The  only  difference  is  that  caution  must  be  used  in 
calculating  the  derivative  matrices,  to  ensure  that  the  appropriate  effective  array  spacing 
is  used  depending  on  whether  (  is  less  than  or  greater  than  L. 

Using  the  results  of  Chapter  5  and  Appendix  A,  the  Cramer-Rao  Bound  on  the 
attitude  estimation  error  using  two  frequencies  is 

-l 


Pcn(q)  = 


2  Re 


(    2L 

£l7d26" 


K)P£{f)6(C) 


l.  C=i 


(7.33) 
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where,  as  with  the  single  frequency, 

4(C)  -  «£)  (m 

d,(C)  =  C^d^C)  (7.35) 

and 

6(C)  =  [MO  MO  d3(C)]  (7.36) 

To  justify  the  observation  made  earlier  that  the  dual-frequency  MLAE  performed  as  well 
or  better  than  the  single  frequency  MLAE  using  data  from  either  frequency,  consider  that 
the  Fisher  information  matrix  for  the  dual-frequency,  FIMDF  is  the  sum  of  FIM  for  each 
frequency.  That  is, 

FIMDF    =    2Re|f;|7c|2D^(C)Pi(g)D(o| 


+        2  Re 


2Re{£|7d2D"(C)Pd1c(9)D(C)j 

£    |7cl26*(0Pi(,,6(0  \  (7-37) 

=    FIMU  +  FIML2  (7.38) 

where  FIML1  and  FIML2  are  the  Fisher  information  matrices  for  the  two  single  frequency 
estimations.  Since  FIMU,  and  FIML2  are  each  always  positive  definite,  FIMDF  must  be 
positive  definite  and  greater2   than  FIMLi  and  FIML2. 

7.5     Summary 
In  this  chapter  we  have  derived  the  dual-frequency  MLAE.  This  estimator  was  shown 
to  optimally  incorporate  information  from  the  two  frequencies  that  carry  the  GPS  P(Y) 
wavefrom.  This  estimator  has  the  same  desirable  properties  of  consistency  and  asymptotic 
efficiency  as  were  demonstrated  for  the  single  frequency  MLAE.  By  examining  the 
Cramer-Rao  Bound  for  the  dual-frequency  MLAE,  it  is  clear  that  it  will  perform  as  well  as 


That  is,  FIMdf  -  FIMU  >  0  and  FIMDF  -  FIML2  >  0. 
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or  better  than  the  single  frequency  MLAE  using  data  from  either  of  the  GPS  frequencies. 
In  addition,  the  dual-frequency  MLAE  significantly  reduces  the  number  of  false  attitudes 
(attitudes  with  metric  values  numerically  close  to  the  value  at  the  true  attitude)  when  the 
sensors  are  spaced  farther  apart  than  the  Nyquist  criterion.  This  gives  it  a  "widelane-like" 
quality,  and  allows  for  better  quality  attitude  estimates  without  introducing  ambiguity. 


CHAPTER  8 
SIMULATIONS  AND  RESULTS 

8.1     Introduction 

This  chapter  presents  quantitative  results  generated  from  a  simulation  of  the  attitude 
estimators  described  in  the  earlier  chapters.  The  components  and  assumptions  of  the 
simulation  are  discussed  below,  followed  by  a  description  and  the  results  from  three 
performance  studies. 

In  this  chapter  abbreviations  for  the  various  attitude  estimators  are  often  used, 
especially  in  the  figures.  Here  we  refer  to  the  Maximum  Likelihood  Attitude  Estima- 
tor algorithm  presented  in  Chapter  4  as  the  "MLAE."  Similarly,  the  direction  finding 
based  approaches  derived  in  Chapter  6  are  referred  to  as  "DFU"  when  the  direction  esti- 
mates1   are  un- weighted  (i.e.  weighting  matrix  is  the  identity  matrix),  and  "DFW"  when 
weighted  by  the  adaptive  SINR.  Finally,  "CONV"  refers  to  attitude  estimation  using  the 
conventional  phase  difference  method  described  in  Chapter  2. 

8.2     Simulation  Methodology 
In  the  context  of  this  work,  a  scenario  constitutes  a  particular  array  topology,  a  set  of 
jammer  locations,  powers,  and  bandwidths,  a  set  of  satellite  locations  and  received  powers, 
the  receiver  integration  time  (time  between  attitude  updates),  and  of  course  the  antenna 
attitude.  The  purpose  of  this  simulation  is  to  demonstrate  the  relative  performance  of 
the  estimators  described  in  the  previous  section  for  various  scenarios.  Specifically,  the 
simulation  is  used  to  evaluate  the  statistical  performance  of  the  estimators  at  a  static 
condition  in  a  Monte  Carlo  fashion.  It  does  not  explore  the  possible  performance  gains 


1  In  the  results  presented  in  this  chapter  both  DFW  and  DFU  use  decoupled  maximum 
likelihood  estimation  to  produce  the  direction  estimates. 
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from  combining  the  GPS  based  attitude  estimates  with  those  of  an  Inertial  Navigation 
System  (INS)  or  other  sensors,  but  focuses  solely  on  the  attitude  derived  from  GPS. 

The  simulation  makes  a  modeling  distinction  between  the  satellite  waveforms  matched 
to  the  despreading  sequence,  which  is  considered  deterministic,  and  the  interference 
which  is  modelled  as  a  random  process.  As  such,  Monte  Carlo  evaluation  of  the  estimator 
performance  for  a  given  scenario  can  be  performed  by  employing  the  estimator  on  the 
superposition  of  the  synthesized  deterministic  data  and  multiple  realizations  of  the 
interference.  Using  this  distinction,  the  simulation  operation  is  as  follows:  first,  the 
"deterministic"  satellite  data  received  at  each  sensor  for  all  satellites  in  view  are  calculated 
using  one  of  the  array  topologies  discussed  below.  Then  the  interference  covariance 
matrices  are  generated,  containing  the  thermal  noise,  contributions  from  each  jammer,  and 
a  model  of  the  multiple  access  interference2   from  other  satellites.  From  the  covariance 
matrix  a  realization  of  the  interference  is  formed  and  added  to  the  satellite  contributions, 
creating  the  demodulated  data  vectors  u(.  An  additional  realization  of  the  interference 
is  created  to  simulate  the  covariance  estimation  process;  each  realization  calculates  its 
own  estimate  of  interference  covariance.  Finally,  the  simulation  implements  the  various 
attitude  estimators  using  the  covariance  matrix  estimates  and  the  received  data  vectors 
calculated  for  the  realization. 

The  relevant  receiver  parameters  are  as  follows.  The  receiver  chain  simulated  has  a 
system  noise  figure  of  4dB,  a  coherent  signal  loss  of  5  dB,  and  a  noncoherent  signal  loss  of 
2  dB.  Each  receiver  hardware  channel  is  assumed  identical.  The  satellite  power  received 
at  the  antenna  is  -163  dB.  The  jammer  effective  radiated  power  (ERP)  is  20  watts,  and 
the  jammers  are  located  20  nautical  miles  from  the  antenna.  Multipath,  either  from  the 
jammer  or  satellites,  is  not  simulated  for  the  studies  below.  As  is  assumed  throughout  this 


2  Since  the  P(Y)  code  multiple  access  interference  is  so  small  compared  to  the  other 
components  of  the  interference,  including  thermal  noise,  in  many  of  the  simulation  studies 
it  is  not  included. 
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work,  all  satellites  used  for  attitude  determination  are  assumed  to  be  in  code  and  Doppler 
track. 

Sensor  Locations 
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Figure  8.1:  Sensor  locations  for  the  three  antenna  topologies.  Quad  (asterisk), Y  (dia- 
mond), and  Hex  (x).  The  Y  antenna  is  actually  a  subset  of  the  Hex  antenna. 


Three  planer  antenna  topologies  are  used  for  the  studies  below.  The  first  antenna 
uses  4  elements  on  a  rectangularly  sampled  grid,  while  the  second  and  third  use  4  and 
7  elements  on  a  triangularly  sampled  grid  in  a  "Y"  and  hexagonal  shape,  respectively. 
Figure  8.1  shows  the  antenna  locations,  in  wavelengths  for  these  antennas.  Except  for 
study  three,  which  uses  wider  baselines,  each  of  these  is  chosen  with  the  maximum  inter- 
element  spacing  allowed  that  prevents  grating  lobes  in  visible  space  [43].  An  interesting 
point  is  that  the  term  "grating  lobes"  used  by  the  antenna  array/radar  community  and 
the  phrase  "integer  ambiguity"  used  in  reference  to  GPS  attitude  determination  are  in 
fact  caused  by  the  same  phenomena.  Both  refer  to  the  situation  where  the  Shannon  (or 
Nyquist)  sampling  criterion  in  the  spatial  domain  is  not  met,  resulting  in  aliasing  (see,  for 
example,  Oppenheim  and  Shafer  [44]). 

Each  sensor  model  incorporates  a  simple  gain  pattern.  Sources  near  boresight  of  the 
sensor  receive  the  most  gain,  while  the  gain  drops  approximately  lOdB  near  the  cutoff 
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point  of  80  degrees  from  boresight.  Satellites  beyond  the  cutoff  point  are  not  included  in 
the  simulation,  since  difficulties  may  exist  with  maintaining  track  on  them.  Figure  8.2 
shows  this  gain  pattern. 


Sensor  Gain  Pattern 
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Figure  8.2:  Sensor  gain  pattern  used  in  the  attitude  simulation. 

For  the  MLAE,  DFW,  and  DFU  simulation  results,  a  modified  alternating  maximiza- 
tion routine  was  employed  to  search  either  the  three  or  two  dimensional  likelihood  space. 
This  approach  was  used  since  computational  speed  was  not  an  issue. 

8.3     Mean  Total  Error 
Each  realization  of  the  simulation  produces  an  attitude  estimate  using  the  estimators 
discussed  in  the  previous  chapters.  There  are  several  ways  in  which  the  simulation  results 
of  this  chapter  may  be  presented.  The  estimated  attitude  could  be  expressed  in  Euler 
angles,  and  statistics  calculated  on  the  errors  in  roll,  pitch,  and  yaw  independently  for 
every  scenario. 

A  more  concise  method,  and  the  one  employed  in  this  work,  is  to  calculate  the  total 
error  for  each  realization.  The  total  error  is  the  angular  rotation  required  to  rotate 
the  antenna  from  the  estimated  attitude  to  the  correct  attitude,  and  is  always  greater 
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than  or  equal  to  zero.  Calculation  of  the  total  error  is  straightforward  using  quaternion 
multiplication.  If  qmeas  represents  the  measured  (estimated)  attitude  quaternion,  and  q 
the  true  attitude  quaternion,  then  the  quaternion  representation  the  rotation  between 
these  two  frames  of  reference,  qerr  is 

qerr  =  q*(qmea*r1  (8.1) 

and  using  the  quaternion  definition  of  equation  (4.34),  where  the  fourth  component  of  the 
quaternion  is  corresponds  to  the  angular  rotation  about  the  axis  defined  by  the  first  three 
components,  the  total  error  is  simply  found  to  be 

Total  Error  =  2  arccos  [qerr(4)]  (8.2) 


For  each  scenario,  the  mean  total  error  across  all  realizations  in  that  scenario  is  calculated 
for  each  of  the  estimation  methods  as  the  average  of  the  total  errors  from  each  simulation 
realization.  Since  the  total  error  is  always  non-negative,  this  method  is  valid  because  one 
error  cannot  algebraically  cancel  another. 

8.4    Study  1:  Single  Jammer  with  "Random"  Location 
In  this  study,  the  estimators  are  evaluated  for  the  three  antenna  topologies  in  a  single 
jammer  scenario.  The  satellite  consteUation  employed  in  the  simulation  uses  the  satellite 
positions  visible  at  latitude  30°  39'  18.70"  North  and  longitude  86°  39'  9.91"  West, 
(located  in  Ft.  Walton  Beach,  Florida)  on  20  April  2001  at  the  time  21:28  UTC.  The 
antenna  orientation  was  in  the  local  level  frame,  resulting  in  7  satellites  being  visible  to 
the  antenna.  In  order  to  evaluate  the  extent  and  variability  of  degradation  possible  from 
a  jammer  at  a  given  displacement  from  boresight,  the  jammer  location  was  moved  across 
the  simulated  sky  with  a  constant  angular  displacement  from  boresight  of  the  array  of 
approximately  44  degrees  (the  projection  onto  the  x-y  plane  of  the  antenna  was  .7).  The 
"randomness"  of  the  jammer  location  at  this  boresight  angle  is  simulated  by  evaluating 
twelve  jammer  positions,  at  0, 30, 60, ... ,  330  degrees  from  North.  Figure  8.3  shows  a 
sine-space  (i.e.  projection  of  the  LOS  vectors  onto  the  plane  of  the  antenna)  plot  of  the 
satellites  and  all  jammer  locations  used  in  this  study. 
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For  each  jammer  location  two  receiver  types  are  evaluated:  one  in  which  the  receiver 
is  configured  to  provide  attitude  estimates  at  50  Hz  and  one  that  provides  estimates  at 
12.5  Hz.  Data  for  the  50  Hz  scenario  are  presented  first,  followed  by  the  data  for  the  12.5 
Hz  scenario.  For  each  data  point,  500  trials  are  simulated. 

12  "Random"  Jammer  Scenarios 


Figure  8.3:  Sine-space  plot  of  the  satellite  (asterisk)  and  the  12  jammer  locations  (X)  for 
the  random  jammer  study.  Only  one  jammer  is  simulated  at  a  time;  in  the  first  scenario 
the  jammer  location  is  at  Ky    =    .7,  and  the  remaining  scenarios  use  the  jammer  positions 
shown  in  order  clockwise.  The  outer  circle  represents  the  horizon,  and  the  area  inside  the 
inner  circle  represents  the  portion  of  the  sky  visible  to  the  sensors. 


It  is  often  useful  to  compare  the  actual  achieved  estimator  performance  to  the 
Cramer-Rao  Bound  (CRB).  For  this  comparison,  we  have  chosen  to  parameterize  attitude 
in  terms  of  Euler  angles,  since  it  seems  more  intuitive  to  discuss  variances  of  Euler  angles 
than  variances  of  quaternion  components.  Figures  8.4,  8.5,  and  8.6  show  the  standard 
deviations  of  MLAE  estimates  in  roll,  pitch,  and  yaw  respectively.  Each  figure  contains 
the  MLAE  estimates  for  the  Quad,  Y,  and  Hex  antennas  as  well  as  the  CRB  for  each 
of  these  antennas  for  the  12  simulated  "random"  jammer  locations.  Notice  that  the 
difference  between  the  MLAE  performance  and  the  CRB  is  in  most  cases  small  and  in 
many  negligible. 
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Standard  Deviation  in  Roll  Channel 
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Figure  8.4:  Standard  deviation  of  antenna  roll  estimates  and  CRB  for  the  Quad  antenna 
(star),  Y  antenna  (square),  and  Hex  antenna  (diamond).  The  12  points  along  the  abscissa 
correspond  to  the  12  jammer  locations  identified  in  Figure  8.3.  The  update  rate  is  50  Hz. 


Standard  Deviation  in  Pitch  Channel 
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Figure  8.5:  Standard  deviation  of  antenna  pitch  estimates  and  CRB  for  the  Quad  antenna 
(star),  Y  antenna  (square),  and  Hex  antenna  (diamond).  The  12  points  along  the  abscissa 
correspond  to  the  12  jammer  locations  identified  in  Figure  8.3.  The  update  rate  is  50  Hz. 
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Standard  Deviation  in  Yaw  Channel 
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Figure  8.6:  Standard  deviation  of  antenna  yaw  estimates  and  CRB  for  the  Quad  antenna 
(star),  Y  antenna  (square),  and  Hex  antenna  (diamond).  The  12  points  along  the  abscissa 
correspond  to  the  12  jammer  locations  identified  in  Figure  8.3.  The  update  rate  is  50  Hz. 
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Insight  into  the  relative  performance  of  the  four  attitude  estimators  discussed  above 
(MLAE,  DFW,  DFU,  and  CONV)  is  provided  by  comparing  the  mean  total  error  for 
each,  for  each  jammer  position.  Figures  8.7,  8.8,  and  8.9  show  the  mean  total  error  these 
estimators  using  the  Quad,  Y,  and  Hex  antennas  respectively. 

Several  observations  may  be  made  about  the  relative  performance  of  these  estimators 
in  the  simulation.  First,  the  Conventional  approach  was  always  the  poorest  performer, 
with  typical  performance  levels  exceeding  8  degrees  of  mean  total  error.  Second,  the 
MLAE  always  provided  the  best  performance.  In  fact,  the  MLAE  in  these  scenarios 
seemed  the  least  sensitive  to  jammer  location  (among  the  12  locations  simulated).  Finally, 
the  DFW  and  DFU  estimators  fell  somewhere  between  the  conventional  and  MLAE. 
Typically,  their  performance  was  much  closer  to  the  MLAE  than  the  conventional,  with 
the  weighting  function  operating  slightly  better  than  the  unweighted. 

Mean  Total  Error  -  1 2  Jam  Locations,  Quad  Ant. 
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Figure  8.7:  Mean  total  angle  error  using  the  Quad  antenna.  The  12  points  along  the 
abscissa  correspond  to  the  12  jammer  locations  identified  in  Figure  8.3.  Performance  is 
shown  for  the  MLAE  (+),  DF-W  (diamond  with  solid  line),  DF-U  (asterisk),  and  conven- 
tional attitude  estimation  algorithms  (diamond  with  dashed  line).  The  update  rate  is  50 
Hz. 
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Mean  Total  Error  -  12  Jam  Locations,  Y  Ant. 
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Figure  8.8:  Mean  total  angle  error  using  the  Y  antenna.  The  12  points  along  the  abscissa 
correspond  to  the  12  jammer  locations  identified  in  Figure  8.3.  Performance  is  shown 
for  the  MLAE  (+),  DF-W  (diamond  with  solid  line),  DF-U  (asterisk),  and  conventional 
attitude  estimation  algorithms  (diamond  with  dashed  line).  The  update  rate  is  50  Hz. 


94 


Mean  Total  Error  -  1 2  Jam  Locations,  Hex  Ant. 
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Figure  8.9:  Mean  total  angle  error  using  the  Hex  antenna.  The  12  points  along  the  ab- 
scissa correspond  to  the  12  jammer  locations  identified  in  Figure  8.3.  Performance  is 
shown  for  the  MLAE  (+),  DF-W  (diamond  with  solid  line),  DF-U  (asterisk),  and  conven- 
tional attitude  estimation  algorithms  (diamond  with  dashed  line).  The  update  rate  is  50 
Hz. 
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Another  way  to  evaluate  these  estimators  is  to  compare  their  performance  in  the 
jammed  environments  to  an  unjammed  case.  Figures  8.10,  8.11,  8.12,  and  8.13  illustrate 
this  comparison  for  the  MLAE,  DFW,  DFU,  and  Conventional  attitude  estimators 
respectively.  Each  plot  presents  data  for  all  three  antenna  topologies. 

An  interesting  observation  is  that  the  performance  of  the  new  estimators  (MLAE, 
DFW,  and  DFU)  in  the  jammed  case  is  not  dramatically  different  than  in  the  non-jammed 
case.  Although  there  is  a  degradation,  the  errors  in  these  single  jammer  scenarios  are 
only  about  twice  as  large  as  those  expected  simply  from  thermal  noise,  and  not  orders  of 
magnitude  larger  as  the  conventional  method  experienced. 

Mean  Total  Error  -  MLAE  with  12  Jam  Locations  and  Unjammed 
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Figure  8.10:  Mean  total  angle  error  using  the  MLAE  algorithm,  comparing  the  single 
jammer  at  "random"  locations  to  the  unjammed  scenario.  The  solid  lines  are  with  the 
jammer,  and  the  dashed  lines  unjammed.  Performance  is  shown  for  the  Hex  antenna  (+), 
Y  antenna  (asterisk),  and  the  Quad  antenna  (diamond).  The  update  rate  is  50  Hz. 
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Mean  Total  Error  -  DF-W  with  12  Jam  Locations  and  Unjammed 
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Figure  8.11:  Mean  total  angle  error  using  the  DF-W  algorithm,  comparing  the  single 
jammer  at  "random"  locations  to  the  unjammed  scenario.  The  solid  lines  are  with  the 
jammer,  and  the  dashed  lines  unjammed.  Performance  is  shown  for  the  Hex  antenna  (+), 
Y  antenna  (asterisk),  and  the  Quad  antenna  (diamond).  The  update  rate  is  50  Hz. 
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Figure  8.12:  Mean  total  angle  error  using  the  DF-U  algorithm,  comparing  the  single 
jammer  at  "random"  locations  to  the  unjammed  scenario.  The  solid  lines  are  with  the 
jammer,  and  the  dashed  lines  unjammed.  Performance  is  shown  for  the  Hex  antenna  (+), 
Y  antenna  (asterisk),  and  the  Quad  antenna  (diamond).  The  update  rate  is  50  Hz. 
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Mean  Total  Error  -  Conv.  with  12  Jam  Locations  and  Unjammed 
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Figure  8.13:  Mean  total  angle  error  using  the  conventional  algorithm,  comparing  the  single 
jammer  at  "random"  locations  to  the  unjammed  scenario.  The  solid  lines  are  with  the 
jammer,  and  the  dashed  lines  unjammed.  Performance  is  shown  for  the  Hex  antenna  (+), 
Y  antenna  (asterisk),  and  the  Quad  antenna  (diamond).  The  update  rate  is  50  Hz. 
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As  the  attitude  update  rate  decreases,  the  amount  of  time  for  satellite  signal  inte- 
gration and  therefore  the  SNR  increases.  The  previous  figures  showed  performance  at  an 
update  rate  of  50  Hz.  For  the  next  set  of  figures  the  update  was  decreased  to  12.5  Hz, 
increasing  integration  time,  SINR,  and  therefore  performance.  Keeping  the  same  format 
as  with  the  50  Hz  data,  figures  8.14,  8.15,  and  8.16  present  the  performance  comparison  to 
the  CRB.  Figures  8.17,  8.18,  and  8.19  present  the  mean  total  error  performance  for  each 
estimator,  and  figures  8.20,  8.20,  8.22,  8.23  contain  the  comparison  between  jammed  and 
unjammed  scenarios. 
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Figure  8.14:  Standard  deviation  of  antenna  roll  estimates  and  CRB  for  the  Quad  antenna 
(star),  Y  antenna  (square),  and  Hex  antenna  (diamond).  The  update  rate  is  12.5  Hz. 
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Standard  Deviation  in  Pitch  Channel 
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Figure  8.15:  Standard  deviation  of  antenna  pitch  estimates  and  CRB  for  the  Quad  an- 
tenna (star),  Y  antenna  (square),  and  Hex  antenna  (diamond).  The  update  rate  is  12.5 
Hz. 
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Figure  8.16:  Standard  deviation  of  antenna  yaw  estimates  and  CRB  for  the  Quad  antenna 
(star),  Y  antenna  (square),  and  Hex  antenna  (diamond).  The  update  rate  is  12.5  Hz. 
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Mean  Total  Error  -  12  Jam  Locations,  Quad  Ant. 
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Figure  8.17:  Mean  total  angle  error  using  the  Quad  antenna.  The  12  points  along  the 
abscissa  correspond  to  the  12  jammer  locations  identified  in  Figure  8.3.  Performance  is 
shown  for  the  MLAE  (+),  DF-W  (diamond  with  solid  line),  DF-U  (asterisk),  and  conven- 
tional attitude  estimation  algorithms  (diamond  with  dashed  line).  The  update  rate  is  12.5 
Hz. 
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Mean  Total  Error  -  1 2  Jam  Locations,  Y  Ant. 


6  8 

Jammer  Location 


10 


12 


Figure  8.18:  Mean  total  angle  error  using  the  Y  antenna.  The  12  points  along  the  abscissa 
correspond  to  the  12  jammer  locations  identified  in  Figure  8.3.  Performance  is  shown 
for  the  MLAE  (+),  DF-W  (diamond  with  solid  line),  DF-U  (asterisk),  and  conventional 
attitude  estimation  algorithms  (diamond  with  dashed  line).  The  update  rate  is  12.5  Hz. 
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Mean  Total  Error  -  1 2  Jam  Locations,  Hex  Ant. 
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Figure  8.19:  Mean  total  angle  error  using  the  Hex  antenna.  The  12  points  along  the  ab- 
scissa correspond  to  the  12  jammer  locations  identified  in  Figure  8.3.  Performance  is 
shown  for  the  MLAE  (+),  DF-W  (diamond  with  solid  line),  DF-U  (asterisk),  and  conven- 
tional attitude  estimation  algorithms  (diamond  with  dashed  line).  The  update  rate  is  12  5 
Hz. 
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Mean  Total  Error  -  MLAE  with  12  Jam  Locations  and  Unjammed 
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Figure  8.20:  Mean  total  angle  error  using  the  MLAE  algorithm,  comparing  the  single 
jammer  at  "random"  locations  to  the  unjammed  scenario.  The  solid  lines  are  with  the 
jammer,  and  the  dashed  lines  unjammed.  Performance  is  shown  for  the  Hex  antenna  (+), 
Y  antenna  (asterisk),  and  the  Quad  antenna  (diamond).  The  update  rate  is  12.5  Hz. 
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Figure  8.21:  Mean  total  angle  error  using  the  DF-W  algorithm,  comparing  the  single 
jammer  at  "random"  locations  to  the  unjammed  scenario.  The  solid  lines  are  with  the 
jammer,  and  the  dashed  lines  unjammed.  Performance  is  shown  for  the  Hex  antenna  (+), 
Y  antenna  (asterisk),  and  the  Quad  antenna  (diamond).  The  update  rate  is  12.5  Hz. 
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Mean  Total  Error  -  DF-U  with  1 2  Jam  Locations  and  Unjammed 
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Figure  8.22:  Mean  total  angle  error  using  the  DF-U  algorithm,  comparing  the  single 
jammer  at  "random"  locations  to  the  unjammed  scenario.  The  solid  lines  are  with  the 
jammer,  and  the  dashed  lines  unjammed.  Performance  is  shown  for  the  Hex  antenna  (+), 
Y  antenna  (asterisk),  and  the  Quad  antenna  (diamond).  The  update  rate  is  12.5  Hz. 
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Figure  8.23:  Mean  total  angle  error  using  the  conventional  algorithm,  comparing  the  single 
jammer  at  "random"  locations  to  the  unjammed  scenario.  The  solid  lines  are  with  the 
jammer,  and  the  dashed  lines  unjammed.  Performance  is  shown  for  the  Hex  antenna  (+), 
Y  antenna  (asterisk),  and  the  Quad  antenna  (diamond).  The  update  rate  is  12.5  Hz. 
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In  the  above  jammed  scenarios,  the  conventional  attitude  estimator  (i.e.  using  phase 
differences  as  the  observables)  was  always  the  worst  performing  algorithm.  As  discussed 
in  Chapter  2,  this  is  because  the  assumptions  in  the  signal  model  are  violated  when  the 
additional  interference  source  is  added  to  the  signals  from  the  satellites.  However,  another 
method  of  comparison  between  these  attitude  estimators  is  in  the  situation  where  the 
assumptions  used  in  the  conventional  model  do  hold,  i.e.  the  unjammed  case. 

Tables  8.1  and  8.2  provide  the  mean  total  error  for  the  unjammed  scenario  for  each 
of  the  three  estimators  developed  in  this  work  and  the  conventional  method.  Now,  unlike 
the  jammed  cases,  the  conventional  approach  outperforms  the  two  approximations  to  the 
MLAE  developed  in  Chapter  6.  This  is  reasonable,  since  the  intermediate  step  of  direction 
finding  involved  in  the  DFW  and  DFU  approaches  may  introduce  error.  However,  the 
conventional  approach  is  still  inferior  to  the  MLAE.  Conceptually,  the  MLAE  outperforms 
the  conventional  approach  by  (optimally)  incorporating  the  amplitude  and  phase  of  the 
response  from  the  satellites,  where  the  conventional  approach3  uses  only  the  phase.  Table 
8.1  presents  the  unjammed  performance  for  the  50  Hz  update  rate,  and  Table  8.2  presents 
the  performance  for  the  12.5  Hz  update  rate.  These  two  rates  correspond  to  the  rates  used 
in  the  previous  plots. 


3  Other  GPS  attitude  determination  works  have  investigated  using  the  amplitude  infor- 
mation as  well  as  the  phase,  but  for  different  reasons  .  In  [1],  the  amplitude  is  used  as  a 
spectral  estimator  and  predictor  of  the  multipath  interference. 
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Table  8.1:  Mean  total  error  performance  comparison  of  the  four  estimators  in  an  un- 
jammed  environment.  Update  rate  is  50  Hz. 


Antenna 

MLAE 

Conventional  AD 

DF  -  Weighted 

DF  -  Unweighted 

Quad 

Y 

Hex 

1.60 
1.13 
0.81 

1.85 
1.26 
0.93 

2.02 
1.50 
1.10 

3.50 
2.79 
2.17 

Table  8.2:  Mean  total  error  performance  comparison  of  the  four  estimators  in  an  un- 
jammed  environment.  Update  rate  is  12.5  Hz. 


Antenna 

MLAE 

Conventional  AD 

DF  -  Weighted 

DF  -  Unweighted 

Quad 

Y 

Hex 

0.78 
0.55 
0.38 

0.92 
0.68 
0.46 

1.07 
0.78 
0.57 

2.18 
1.62 
1.18 
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8.5     Study  2:  Varying  Attitude  Update  Rate 
In  attitude  determination,  the  attitude  should  be  essentially  constant  during  an 
integration  period  to  prevent  smearing  of  attitude  data.  Also,  the  SNR  (and  measurement 
accuracy)  increases  with  integration  time.  Therefore  the  update  rate  is  a  key  design 
parameter  in  attitude  determination,  and  involves  engineering  tradeoffs.  More  frequent 
estimates  serve  to  ensure  that  the  attitude  is  indeed  constant  during  the  estimation 
process,  but  have  a  larger  variance  since  less  time  is  available  for  integration  of  the 
GPS  signals.  Similarly,  less  frequent  updates  are  more  accurate,  but  only  if  the  antenna 
attitude  is  essentially  constant  during  the  entire  integration  process.  Additional  antennas 
can  be  used  to  improve  the  accuracy,  but  at  the  cost  of  increased  system  complexity.  The 
update  rate  is  therefore  a  compromise  between  expected  platform  dynamics,  accuracy,  and 
complexity. 

In  this  study,  we  examine  the  performance  as  a  function  of  integration  time.  Case  1 
evaluates  this  performance  for  a  single  jammer,  and  Case  2  for  3  jammers  in  the  field  of 
view.  As  in  Study  1,  the  jammers  are  located  at  an  angle  of  44  -  48  degrees  from  antenna 
boresight,  as  shown  in  figure  8.24.  In  this  scenario,  the  antenna  attitude  is  non-zero  only 
in  yaw,  with  £*  =  § . 

As  with  the  random  jammer  study,  MLAE  performance  is  compared  to  the  CRB. 
Figures  8.25,  8.26,  and  8.27  present  the  standard  deviations  of  MLAE  estimates  in  roll, 
pitch,  and  yaw  respectively,  as  a  function  of  update  rate,  for  all  three  antenna  topologies. 
Recall  from  Chapter  5  that  the  MLAE  is  asymptotically  efficient.  This  is  equivalent  to 
stating  that  as  the  effective  SNR  increases,  the  MLAE  approaches  the  CRB.  There  are 
two  ways  in  which  the  effective  SNR  may  increase:  as  the  integration  time  increases 
(update  rate  decreases),  and  as  the  number  and  /  or  spacing  of  sensors  increases.  The 
Hex  antenna  incorporates  more  sensors  than  the  Quad  or  Y  antenna,  and  the  difference 
between  the  MLAE  standard  deviations  and  the  CRB  with  this  antenna  are  small 
throughout  the  span  of  update  rates  investigated.  The  antenna  with  the  smallest  spacing 
and  number  of  sensors  is  the  Quad.  For  fast  update  rates  a  discernable  difference  between 
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the  MLAE  performance  and  the  CRB  exists  with  this  antenna  topology,  however  this 
difference  decreases  as  the  integration  time  increases. 


Jammer  /  Satellite  Locations 


Figure  8.24:  Satellite  (asterisk)  and  Jammer  Locations  (x)  for  the  varying  update  rate 
study.  The  jammer  indicated  by  the  arrow  is  the  only  jammer  used  in  case  1.  In  case  2, 
all  three  jammers  appear. 
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Figure  8.25:  Standard  deviation  of  antenna  roll  estimates  and  CRB  for  the  Quad  antenna 
(star),  Y  antenna  (square),  and  Hex  antenna  (diamond)  vs.  update  rate.  One  jammer  in 
view. 
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Figure  8.26:  Standard  deviation  of  antenna  pitch  estimates  and  CRB  for  the  Quad  an- 
tenna (star),  Y  antenna  (square),  and  Hex  antenna  (diamond)  vs.  update  rate.  One 
jammer  in  view. 
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Standard  Deviation  in  Yaw  Channel 
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Figure  8.27:  Standard  deviation  of  antenna  yaw  estimates  and  CRB  for  the  Quad  antenna 
(star),  Y  antenna  (square),  and  Hex  antenna  (diamond)  vs.  update  rate.  One  jammer  in 
view. 
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It  is  insightful  to  compare  the  relative  performance  of  the  four  attitude  estimators  by 
examining  their  mean  total  error  as  the  update  rate  changes  for  a  fixed  jammer  scenario. 
Figures  8.28,  8.29,  and  8.30  present  this  comparison  for  the  Quad,  Y,  and  Hex  antennas 
respectively. 

The  performance  for  the  Hex  antenna  topology,  which  has  7  sensors,  is  superior 
to  the  4  sensor  Y  and  Quad  topologies.  The  Y,  which  is  set  on  the  triangular  lattice, 
outperforms  the  rectangular  lattice  Quad  antenna  because  the  sensors  are  slightly  further 
apart  for  the  Y,  resulting  in  longer  baseline  lengths.  The  MLAE  estimator  always 
outperforms  all  other  attitude  estimators.  Performance  for  any  estimator  or  antenna 
increases  as  the  update  rate  decreases,  since  the  resulting  longer  integration  times  increase 
the  signal  to  interference  ratio.  As  the  integration  increases,  the  performance  of  even 
conventional  attitude  determination  methods  improves,  however  its  error  is  several  times 
greater  than  the  methods  developed  in  this  work  even  at  the  relative  slow  attitude  update 
rates  of  10  Hz. 
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Figure  8.28:  Mean  total  angle  error  using  the  Quad  antenna  vs.  update  rate  for  1  jam- 
mer. Performance  is  shown  for  the  MLAE  (+),  DF-W  (diamond  with  solid  line),  DF-U 
(asterisk),  and  conventional  attitude  estimation  algorithms  (diamond  with  dashed  line). 
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Mean  Total  Error,  1  Jammer,  Y  Ant. 
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Figure  8.29:  Mean  total  angle  error  using  the  Y  antenna  vs.  update  rate  for  1  jammer. 
Performance  is  shown  for  the  MLAE  (+),  DF-W  (diamond  with  solid  line),  DF-U  (aster- 
isk), and  conventional  attitude  estimation  algorithms  (diamond  with  dashed  line). 
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Mean  Total  Error,  1  Jammer,  Hex  Ant. 
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Figure  8.30:  Mean  total  angle  error  using  the  Hex  antenna  vs.  update  rate  for  1  jam- 
mer. Performance  is  shown  for  the  MLAE  (+),  DF-W  (diamond  with  solid  line),  DF-U 
(asterisk),  and  conventional  attitude  estimation  algorithms  (diamond  with  dashed  line). 
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Figures  8.31,  8.32,  8.33  and  8.34  compare  the  performance  of  the  MLAE,  DFW, 

DFU  and  conventional  attitude  estimators  in  the  single  jammer  environment  to  their 

performance  in  an  unjammed  environment.  This  provides  an  indication  as  to  the  level  of 

robustness,  or  jammer  tolerance  of  each  algorithm.  The  MLAE  has  the  greatest  jammer 

tolerance.  The  maximum  increase  in  mean  total  error  between  jammed  and  unjammed 

is  a  factor  of  around  1.5  for  the  Quad  antenna  and  much  less  for  the  other  two.  The 

DFW  maximum  increase  is  greater  than  three  for  the  Quad  antenna,  and  less  for  the 

others  as  well.  It  is  interesting  that  the  DFU  maximum  increase  is  less  than  the  DFW, 

but  in  general  its  ratio  of  jammed  to  unjammed  mean  total  error  is  larger  than  that  of 

the  DFW,  and  the  unjammed  performance  of  the  DFU  is  worse  than  the  DFW.  Finally, 

the  conventional  approach  has  the  largest  increase  in  mean  total  error  in  a  jammed 

environment  over  an  unjammed  scenario.  Of  course,  this  was  expected  and  was  in  fact  the 

motivation  behind  this  research. 
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Figure  8.31:  Mean  total  angle  error  using  the  MLAE  algorithm,  comparing  the  single  jam- 
mer at  "random"  locations  to  the  unjammed  scenario,  vs.  update  rate.  The  solid  lines 
are  with  the  jammer,  and  the  dashed  lines  unjammed.  Performance  is  shown  for  the  Hex 
antenna  (+),  Y  antenna  (asterisk),  and  the  Quad  antenna  (diamond). 
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Mean  Total  Error  -  DF-W  with  1  Jammer  and  Unjammed 
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Figure  8.32:  Mean  total  angle  error  using  the  DF-W  algorithm,  comparing  the  single 
jammer  scenario  to  the  unjammed  scenario,  vs.  update  rate.  The  solid  lines  are  with  the 
jammer,  and  the  dashed  lines  unjammed.  Performance  is  shown  for  the  Hex  antenna  (+), 
Y  antenna  (asterisk),  and  the  Quad  antenna  (diamond). 
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Mean  Total  Error  -  DF-U  with  1  Jammer  and  Unjammed 
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Figure  8.33:  Mean  total  angle  error  using  the  DF-U  algorithm,  comparing  the  single  jam- 
mer scenario  to  the  unjammed  scenario,  vs.  update  rate.  The  solid  lines  are  with  the 
jammer,  and  the  dashed  lines  unjammed.  Performance  is  shown  for  the  Hex  antenna  (+), 
Y  antenna  (asterisk),  and  the  Quad  antenna  (diamond). 
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Mean  Total  Error  -  Conv.  with  1  Jammer  and  Unjammed 
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Figure  8.34:  Mean  total  angle  error  using  the  conventional  algorithm,  comparing  the  single 
jammer  scenario  to  the  unjammed  scenario,  vs.  update  rate.  The  solid  lines  are  with  the 
jammer,  and  the  dashed  lines  unjammed.  Performance  is  shown  for  the  Hex  antenna  (+), 
Y  antenna  (asterisk),  and  the  Quad  antenna  (diamond). 
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Figures  8.35,  8.36,  and  8.37  present  the  mean  total  error  for  the  four  attitude  estima- 
tors in  a  three  jammer  scenario.  This  scenario  with  three  jammers  is  stressing  for  the  four 
sensor  antenna  topologies,  and  this  is  reflected  in  the  performance  estimates.  Not  every 
scenario  involving  three  jammers  will  exhibit  the  performance  reported  here.  In  fact,  the 
three  jammers  are  spaced  throughout  the  sky,  in  essentially  a  "worst  case"  type  scenario. 
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Figure  8.35:  Mean  total  angle  error  using  the  Quad  antenna  vs.  update  rate  for  3  jam- 
mers. Performance  is  shown  for  the  MLAE  (+),  DF-W  (diamond  with  solid  line),  DF-U 
(asterisk),  and  conventional  attitude  estimation  algorithms  (diamond  with  dashed  line). 
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Mean  Total  Error,  3  Jammers,  Y  Ant. 
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Figure  8.36:  Mean  total  angle  error  using  the  Y  antenna  vs.  update  rate  for  3  jammers. 
Performance  is  shown  for  the  MLAE  (+),  DF-W  (diamond  with  solid  line),  DF-U  (aster- 
isk), and  conventional  attitude  estimation  algorithms  (diamond  with  dashed  line). 
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Mean  Total  Error,  3  Jammers,  Hex  Ant. 
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Figure  8.37:  Mean  total  angle  error  using  the  Hex  antenna  vs.  update  rate  for  3  jam- 
mers. Performance  is  shown  for  the  MLAE  (+),  DF-W  (diamond  with  solid  line),  DF-U 
(asterisk),  and  conventional  attitude  estimation  algorithms  (diamond  with  dashed  line). 
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Figures  8.38,  8.39,  8.40,  and  8.41  present  the  performance  of  the  MLAE,  DFW, 
DFU  and  conventional  attitude  estimators  in  the  three  jammer  environment  and  their 
performance  in  an  unjammed  environment.  In  all  cases  the  ratio  between  the  jammed 
and  unjammed  cases  is  larger  in  this  three  jammer  case  than  in  that  of  the  single  jammer. 
However,  the  same  relative  robustness  exists  among  the  four  estimators,  with  the  MLAE 
being  the  most  robust  and  the  conventional  approach  the  least. 
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Figure  8.38:  Mean  total  angle  error  using  the  MLAE  algorithm,  comparing  the  three 
jammer  scenario  to  the  unjammed  scenario,  vs.  update  rate.  The  solid  lines  are  with  the 
jammer,  and  the  dashed  lines  unjammed.  Performance  is  shown  for  the  Hex  antenna  (+), 
Y  antenna  (asterisk),  and  the  Quad  antenna  (diamond). 
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Figure  8.39:  Mean  total  angle  error  using  the  DF-W  algorithm,  comparing  the  three  jam- 
mer scenario  to  the  unjammed  scenario,  vs.  update  rate.  The  solid  lines  are  with  the 
jammer,  and  the  dashed  lines  unjammed.  Performance  is  shown  for  the  Hex  antenna  (+), 
Y  antenna  (asterisk),  and  the  Quad  antenna  (diamond). 
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Figure  8.40:  Mean  total  angle  error  using  the  DF-U  algorithm,  comparing  the  three  jam- 
mer scenario  to  the  unjammed  scenario,  vs.  update  rate.  The  solid  lines  are  with  the 
jammer,  and  the  dashed  lines  unjammed.  Performance  is  shown  for  the  Hex  antenna  (+), 
Y  antenna  (asterisk),  and  the  Quad  antenna  (diamond). 
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Mean  Total  Error  -  Conv.  with  3  Jammers  and  Unjammed 
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Figure  8.41:  Mean  total  angle  error  using  the  conventional  algorithm,  comparing  the  three 
jammer  scenario  to  the  unjammed  scenario,  vs.  update  rate.  The  solid  lines  are  with  the 
jammer,  and  the  dashed  lines  unjammed.  Performance  is  shown  for  the  Hex  antenna  (+), 
Y  antenna  (asterisk),  and  the  Quad  antenna  (diamond). 
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8.6     Study  3:  Wider  Baselines  and  Dual  Frequency 


In  this  study  we  will  show  performance  results  from  simulations  of  the  dual-frequency 
attitude  estimator.  In  addition,  this  study  considers  a  larger  baseline  length  than  the 
Nyquist  spacing  used  in  the  previous  studies.  The  sensors  are  spaced  at  five  times  the 
spacing  used  before,  which  corresponds  to  approximately  .6  meters  apart  for  the  quad 
antenna,  and  slightly  larger  for  the  other  two  antennas.  Figure  8.42  shows  the  sensor 
locations  used  in  this  study.  This  is  why  the  performance  is  so  dramatically  better  than  in 
the  two  previous  studies. 

This  study  does  not  address  the  decrease  in  false  attitudes  discussed  in  Chapter 
7,  but  only  the  increase  in  performance  attained  by  including  the  second  frequency.  In 
addition,  like  all  the  previous  studies,  assumed  that  the  adaptive  antenna  array  is  able  to 
keep  the  satellites  in  code  and  Doppler  track  even  in  the  presence  of  jammers. 

Figures  8.43,  8.44,  8.45  present  the  mean  total  error  performance  of  the  single 
frequency  MLAE  operating  at  the  LI  GPS  frequency,  the  L2  GPS  frequency,  and  the 
dual-frequency  MLAE  which  uses  both  frequencies  .  The  scenario  is  the  single  jammer 
scenario  used  in  Study  2  Case  1.  Figures  8.46,  8.47,  8.48  present  similar  plots  but  using 
the  three  jammer  scenario  of  Study  2  Case  2.  However  in  the  scenario  that  provided  the 
data  for  all  of  these  figures  the  jammers  were  present  on  both  GPS  frequencies  at  the  same 
power  levels  per  band  as  used  before. 

Some  observations  may  be  drawn  from  these  simulation  results.  When  the  perfor- 
mance is  noticeably  better  at  one  frequency  than  the  other  (due  to  the  complex  grating 
lobe  structure  and  the  jammer  location),  the  dual-frequency  MLAE  performs  slightly 
better  than  the  best  performing  single  frequency.  When  the  two  single  frequency  MLAE's 
perform  about  the  same,  the  dual-frequency  MLAE  performs  noticeably  better  than  the 
two  single  frequencies.  In  all  cases  the  dual-frequency  MLAE  performed  as  well  as  or 
better  than  the  single  frequency  MLAE  operating  on  either  band.  These  observations  are 
consistent  with  the  performance  improvements  and  CRB  derivation  discussed  in  Chapter 
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Figure  8.42:  Sensor  Locations  used  in  the  wide  baseline  study.  The  three  antenna  topolo- 
gies are  the  Quad  (asterisk), Y  (diamond),  and  Hex  (x). 
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Figure  8.43:  Mean  total  angle  error  using  the  Quad  antenna  vs.  update  rate  for  1  jammer. 
Performance  is  shown  for  the  single  frequency  MLAE  at  LI  (diamond)  and  L2  (+),  and 
for  the  dual  frequency  MLAE  (asterisk). 
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Figure  8.44:  Mean  total  angle  error  using  the  Y  antenna  vs.  update  rate  for  1  jammer. 
Performance  is  shown  for  the  single  frequency  MLAE  at  LI  (diamond)  and  L2  (+),  and 
for  the  dual  frequency  MLAE  (asterisk). 
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Figure  8.45:  Mean  total  angle  error  using  the  Hex  antenna  vs.  update  rate  for  1  jammer. 
Performance  is  shown  for  the  single  frequency  MLAE  at  Ll  (diamond)  and  L2  (+),  and 
for  the  dual  frequency  MLAE  (asterisk). 
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Mean  Total  Error,  3  Jammers,  Wide  Baseline,  Quad  Ant. 
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Figure  8.46:  Mean  total  angle  error  using  the  Quad  antenna  vs.  update  rate  for  3  jam- 
mers. Performance  is  shown  for  the  single  frequency  MLAE  at  Ll  (diamond)  and  L2  (+), 
and  for  the  dual  frequency  MLAE  (asterisk). 
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Figure  8.47:  Mean  total  angle  error  using  the  Y  antenna  vs.  update  rate  for  3  jammers. 
Performance  is  shown  for  the  single  frequency  MLAE  at  Ll  (diamond)  and  L2  (+),  and 
for  the  dual  frequency  MLAE  (asterisk). 
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Figure  8.48:  Mean  total  angle  error  using  the  Hex  antenna  vs.  update  rate  for  3  jammers. 
Performance  is  shown  for  the  single  frequency  MLAE  at  Ll  (diamond)  and  L2  (+),  and 
for  the  dual  frequency  MLAE  (asterisk). 
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8.7     Conclusions 


In  this  chapter  we  have  demonstrated  via  simulation  the  performance  of  the  three 
new  attitude  estimators  developed  in  the  previous  chapters.  Three  separate  studies 
were  undertaken,  and  performance  of  the  estimators  was  compared  to  each  other,  the 
conventional  attitude  determination  method,  the  Cramer-Rao  Bound,  and  to  unjammed 
cases  for  three  antenna  topologies. 

Several  observations  were  made  from  these  studies  that  are  summarized  here.  First, 
the  MLAE  always  provides  the  best  mean  total  attitude  error  performance.  This  is 
true  regardless  of  the  antenna  topology,  number  of  jammers,  or  update  rate.  Second, 
the  MLAE  was  shown  to  approach  (or  essentially  equal)  the  CRB.  Since  the  MLAE  is 
asymptotically  efficient,  the  difference  between  MLAE  performance  and  the  CRB  depends 
on  the  number  of  sensors,  number  of  jammers,  array  topology,  and  update  rate.  In 
general,  the  stronger  the  satellite  signal  in  comparison  to  the  interference,  the  smaller  the 
difference  between  the  MLAE  and  the  CRB. 

Third,  the  two  approximations  to  the  MLAE,  i.e.  the  DFW  and  DFU  algorithms, 
often  estimate  attitude  with  errors  close  to  but  slightly  larger  than  the  MLAE,  with  DFW 
typically  outperforming  DFU.  This  is  important,  since  the  computational  savings  from 
these  suboptimal  algorithms  may  be  significant  in  a  tactical  application. 

Fourth,  the  conventional  attitude  estimation  approach  (using  phase  differences)  has 
very  little  performance  in  a  jammed  environment.  As  discussed  earlier,  this  is  expected 
and  provided  the  motivation  to  investigate  jam-resistant  algorithms.  It  is  interesting  to 
note,  however,  that  in  an  unjammed  environment  the  conventional  approach  outperforms 
all  other  estimators  except  the  MLAE.  This  implies  that  a  tactical  approach  using  the 
DFW  or  DFU  methods  might  need  to  investigate  switching  between  the  direction  finding 
algorithms  and  the  conventional  approach  depending  on  the  presence  or  absence  of 
interference  sources. 

Finally,  the  study  involving  the  dual-frequency  MLAE  showed  that  it  always  outper- 
formed the  single  frequency  MLAE  at  either  band.  The  amount  of  the  performance  gain  is 
driven  by  the  relative  performance  in  each  band.  That  is,  if  one  band  is  providing  poorer 
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performance  that  the  other,  the  performance  gain  of  the  dual-frequency  MLAE  over  the 
single  frequency  MLAE  at  the  better  performing  band  will  be  small. 


CHAPTER  9 
CONCLUSIONS 

9.1     Summary 

This  dissertation  provides  a  comprehensive  investigation  into  the  new  field  of  robust, 
jam-resistant  attitude  determination  using  the  Global  Positioning  System.  With  its  low 
transmit  power  and  considerable  distance  from  receivers,  GPS  is  known  to  be  susceptible 
to  jamming.  Efforts  to  date  have  investigated  various  receiver  designs  and  adaptive 
antenna  arrays  to  improve  the  anti-jam  capabilities  of  GPS  position  location,  but  until  this 
work  no  extension  of  anti-jam  capabilities  to  attitude  determination  existed. 

Direction  finding  and  GPS  based  attitude  determination  are  similar  concepts  related 
by  their  use  of  carrier  phase  interferometry.  Historically,  direction  finding  algorithms  cast 
in  terms  of  a  maximum  likelihood  estimator  have  shown  significant  interference  resistance, 
so  the  approach  of  this  work  was  to  exploit  similarities  between  these  estimators  and  atti- 
tude determination.  Recent  works  in  multiple  source  direction  finding  have  demonstrated 
that  if  the  source  waveforms  are  uncorrelated,  the  maximum  likelihood  estimates  of  all 
directions  are  (asymptotically)  equal  to  decoupled  estimates  of  individual  directions,  i.e. 
the  direction  to  one  source  may  be  found  without  knowledge  or  use  of  the  direction  to 
any  other.  This  is  an  important  result  for  the  direction  finding  application,  and  provided 
key  insight  for  this  dissertation.  However,  use  of  the  sources  independently  does  not  lead 
to  the  best  possible  performance  for  the  attitude  determination  application.  Since  the 
positions  of  the  sources  (the  GPS  satellites)  with  respect  to  each  other  are  known,  to 
treat  each  source  independently  and  to  not  make  use  of  this  information  would  lead  to  a 
sub-optimal  solution. 

An  optimal  solution  does  exist  if  some  of  the  concepts  from  the  direction  finding  field 
are  redefined.  Specifically,  by  reparameterizing  the  argument  of  the  array  response  vector 
to  be  antenna  attitude  and  local  level  frame  direction  instead  of  direction  in  the  antenna 
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frame,  a  coupled  estimator  is  derived  that  makes  use  of  the  known  relative  positions  of  the 
satellites.  This  approach  is  named  the  Maximum  Likelihood  Attitude  Estimator  (MLAE). 
By  using  estimates  of  the  interference  statistics,  this  estimator  reduces  to  a  simple  search 
over  the  uncertainty  in  (the  three  parameters  of)  attitude  of  the  maximum  likelihood 
metric.  The  estimator  decomposes  into  a  sum  of  terms,  where  each  term  is  composed  of 
data  relating  to  a  single  satellite,  providing  for  a  parallel  implementation  architecture. 

In  this  dissertation  the  MLAE  is  shown  to  possess  several  desirable  statistical  proper- 
ties. The  MLAE  is  proven  to  be  a  consistent  estimator.  That  is,  as  the  integration  time 
increases  the  variance  of  the  estimation  error  decreases.  In  addition,  it  is  asymptotically 
statistically  efficient.  As  the  integration  time  increases,  the  estimation  error  approaches 
the  lower  limit  of  any  unbiased  estimator. 

In  addition  to  the  MLAE,  this  work  considers  two  approximations  to  this  algorithm. 
These  approaches  reduce  the  computational  burden  of  searching  a  three  dimensional 
attitude  space  by  finding  the  directions  to  the  satellites  independently,  which  is  typically 
a  series  of  independent  two  dimensional  searches,  one  per  source.  The  difference  between 
the  two  approximations  is  in  the  conversion  from  the  direction  estimates  back  to  the 
attitude.  In  the  first  algorithm,  the  contribution  from  each  satellite  is  weighted  equally, 
while  in  the  second  each  satellite  is  weighted  by  its  adaptive  signal  to  interference  ratio 
(SINR).  Simulation  results  indicate  that,  in  general,  inclusion  of  the  adaptive  SINR 
weighting  improves  performance  of  the  estimator  with  a  very  small  increase  in  the  number 
of  computations. 

Currently  each  satellite  broadcasts  the  GPS  P(Y)  code  on  two  frequencies.  This 
provides  additional  information  that  can  be  incorporated  to  improve  the  quality  of  the 
attitude  estimate.  To  exploit  this  information,  the  dissertation  develops  the  dual-frequency 
MLAE.  The  dual-frequency  MLAE  is  also  asymptotically  statistically  efficient,  and 
provides  two  benefits  over  the  single  frequency  version.  First,  the  estimation  error  is 
reduced  by  incorporating  the  second  frequency.  Second,  when  only  a  few  satellites  are 
being  tracked  and  used  for  attitude  estimation  and  the  sensor  spacing  is  greater  than  the 
Nyquist  criterion,  incorporating  the  second  frequency  "smooths"  the  likelihood  surface, 
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providing  less  areas  that  may  lead  to  incorrect  solutions.  This  important  result  provides 
robustness  against  false  solutions  when  the  sensor  spacing  is  increased,  which  is  known  to 
increase  estimator  accuracy. 

A  series  of  performance  studies  were  performed  and  documented  in  this  dissertation. 
These  studies  compared  the  performance  of  the  various  attitude  estimators  in  a  variety 
of  jammed  environments  to  one  another,  to  the  Cramer-Rao  Bound  (CRB),  and  to  per- 
formance in  an  unjammed  case.  In  all  cases  the  MLAE  provided  superior  performance  to 
all  other  attitude  estimators  simulated.  Furthermore,  the  MLAE  and  its  approximations 
significantly  outperformed  the  conventional  approach  in  all  jammed  cases  investigated. 

9.2     Future  Work 

This  research  has  taken  important  steps  towards  understanding  and  mitigating 
jamming  and  other  interference  for  GPS  attitude  determination.  There  are,  however, 
additional  contributions  beyond  this  research  that  wiU  increase  the  performance  of  these 
methods,  provide  better  estimates  of  sensitivity  to  real  world  limitations,  or  both.  A  few 
of  these  contributions  are  listed  below. 

As  discussed  in  Chapter  4,  the  MLAE  involves  a  three  dimensional  search  over 
attitude.  This  can  be  computationally  expensive,  and  may  be  prohibitive  for  large 
uncertainties  in  attitude.  Although  several  maximization/minimization  techniques 
exist,  a  comparison  or  optimization  of  them  has  not  been  performed  for  this  particular 
application.  In  addition,  it  may  be  possible  to  develop  a  reasonably  accurate  model  for 
the  likelihood  volume,  and  with  a  few  hundred  data  points  estimate  the  parameters  of  this 
model.  If  the  parameters  are  estimated  with  sufficient  accuracy,  the  model  will  provide  a 
good  estimate  of  the  location  of  the  extrema.  From  this  model  estimate,  either  a  fine  grid 
search  or  simply  a  check  of  neighboring  points  could  be  performed  to  produce  the  final 
attitude  estimate. 

Several  assumptions  which  tend  to  idealize  performance  were  made  concerning  the 
antenna  and  receiver  characteristics.  Mutual  coupling  between  sensors  was  ignored, 
and  the  receiver  channels  were  assumed  to  be  perfectly  matched.  These  "real  world" 
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imperfections  will  tend  to  degrade  performance.  However,  this  effect  can  be  mitigated  to 
some  degree  by  levying  appropriate  specifications  on  the  hardware  design. 

Another  real  world  area  of  investigation  involves  error  in  the  array  calibration.  The 
simulation  results  provided  in  this  dissertation  assume  that  the  array  is  calibrated  (i.e. 
the  array  response  vector  is  known  for  all  spatial  locations  of  interest).  However,  if  the 
effective  locations  (phase  centers)  of  the  sensors  are  in  error,  the  array  response  vector 
will  not  be  correct.  This  will  degrade  performance  for  both  the  MLAE  and  the  direction 
finding  based  approaches.  One  investigation  would  be  to  study  the  relative  robustness 
of  the  MLAE  and  the  direction  finding  attitude  determination  approaches  to  array 
calibration  errors. 

In  this  research,  simulation  results  have  focused  on  the  statistics  (either  standard 
deviation  in  Euler  angle  measurements  or  mean  total  error)  of  the  estimators  at  a  static 
snapshot  in  time.  In  addition,  an  assumption  on  results  has  been  that  the  adaptive 
beamformer  portion  (the  typical  GPS  anti-jam  system  for  position  location)  has  converged 
to  keep  the  satellites  in  track.  The  next  level  of  investigation  of  these  concepts  is  to 
explore  the  performance  of  a  closed  loop  system.  This  could  include  a  tracker  design  that 
incorporates  both  the  GPS  based  attitude  information  and  an  inertial  navigation  system 
(INS)  that  uses  gyros  for  attitude  estimation.  Inherent  in  this  design  would  be  a  method 
to  estimate  the  attitude  errors  from  the  new  algorithms,  possibly  via  the  Cramer  Rao 
Bound.  This  investigation  could  also  include  antenna  (body)  motion.  This  may  have  the 
effect  of  smearing  data  during  the  observation  interval,  and  would  allow  the  beamformer 
to  adaptively  attempt  to  remove  jammers  to  keep  the  satellites  in  track  in  code  and 
Doppler. 

The  derivations  in  this  work  do  not  include  multipath  reflections  of  either  the  jammer 
waveforms  or  the  satellites.  To  account  for  these,  multiple  snapshots  of  demodulated 
data  may  be  used  to  optimally  remove  interference  and  estimate  the  attitude.  This 
would  essentially  be  a  Space  -  Time  Adaptive  Processing  (STAP)  extension  to  attitude 
determination.  Research  in  this  area  may  focus  on  quantifying  the  extent  of  possible 
performance  gains,  and  the  best  integration  of  STAP  and  code  tracking  loops. 
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Finally,  further  research  is  warranted  in  the  dual  frequency  work  Chapter  7.  The  dual 
frequency  extension  of  the  MLAE  developed  in  that  chapter  may  be  just  one  example  in  a 
class  of  multi-spectral  estimation  algorithms. 


APPENDIX  A 
CRAMEr-RAO  BOUND 

The  Cramer-Rao  Bound  (CRB)  provides  a  lower  bound  on  the  mean  squared  error 
(MSE)  of  the  estimation  error.  An  estimator  whose  MSE  achieves  the  CRB  is  said  to  be 
statistically  efficient  (see,  for  example,  Stoica  and  Moses  [9]).  In  this  section  we  will  derive 
the  CRB  for  the  attitude  estimation  problem. 

The  Cramer-Rao  Bound  (CRB)  matrix  P„.  is  found  from  [9]  as 


P„=     E 


«n/x(xk>) 


d\nfx(X\P) 


dp 


(A.l) 


dp 

where  p  is  a  vector  of  the  real  unknown  parameters  (  i.e.  the  parameters  to  be  estimated), 
and  fx(x\p)  is  the  likelihood  function.  Before  performing  the  matrix  inversion,  the  right 
side  of  the  equation,  before  performing  the  inversion  (i.e.  P"1),  is  the  "Fisher  Information 
Matrix,"  and  will  be  denoted  here  by  FIM.  For  the  attitude  determination  when  the 
interference  covariance  is  known  and  each  satellite  return  is  scaled  by  an  unknown 
complex  gain,  p  consists  of  the  three  unknown  parameters  of  attitude,  qu  i  =  1, 2, 3  and 
the  2  x  L  real  components  of  the  L  complex  gains.  It  is  useful  to  organize  these  terms  in 
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the  following  order: 


P  = 


(A.2) 


<h 

Re(7i) 
Re(72) 

Re(7L) 
Im(7i) 
Im(72) 


Im(7t) 

A  convenient  parameterization  of  the  three  attitude  parameters  for  the  CRB  calculation, 
and  the  one  used  in  this  work,  is  to  choose  the  three  Euler  angles  £$,  £e,  and  £*. 
Using  (4.49),  the  log-likelihood  ratio  (LLR)  may  be  written  as 

L 

lnA(x|p)  =  {constant  terms}  -  J^[u(  -  ^ajfa^Cf^ui  -  -yiau(q)}  (A.3) 

i=i 

In  order  to  evaluate  (A.l)  (i.e.  to  generate  the  various  terms  of  the  FIM),  the  following 

three  derivative  formulas  are  required: 


«n/x(x|p) 
0Re(7«) 


2Re{a?((7)Cr1wi} 


(A.4) 


<9Im(7() 


=  2Im{af(g)Cr1wi} 


dlnfx(X\p) 
dqi 


=  2  Re 


S^df^crvj 


(A.5) 


(A.6) 


where  w(  is  defined  in  (4.35)  and  dk(l)  is  defined  as  the  derivative  of  a;  with  respect  to  the 

kth.  term  of  attitude  q,  qk. 

.    ...        dan 

(A.7) 


oqk 


— 
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Note  that  each  of  the  L  array  response  vectors  corresponding  to  the  L  satellites  will  have 
different  derivative  vectors.  The  convention  used  here  is  as  follows:  the  subscript  indicates 
the  attitude  parameter  with  respect  to  which  the  derivative  is  being  taken,  while  the  value 
in  parentheses  indicates  the  satellite  to  which  the  array  response  vector  corresponds.  For 
notational  convenience,  a  whitened  form  of  the  derivative  vectors  is  defined  as 

dt(0  =  Cr1/2d;  (A.8) 

where 

cr1  =  c;1/2c;1/2  (A.9) 


.-1/2 


=  (cr1/2)"  (a.io) 


The  order  of  the  terms  in  the  vector  p  determines  the  structure  of  the  FIM.  It  is  easy 
to  see  that  p  is  composed  of  three  types  of  terms:  "attitude"  terms,  "real  gain"  terms, 
and  "imaginary  gain"  terms.  The  FIM  will  then  be  composed  of  blocks  of  (the  expected 
value  of)  the  products  of  the  derivatives  with  respect  to  each  of  the  types  of  terms.  For 
example,  the  upper  3x3  block  of  the  FIM  is  the  "attitude-attitude"  block,  and  the 
3  x  L  block  immediately  to  the  right  is  the  "attitude-real  gain"  block,  and  so  on.  The 
approach  employed  here  for  populating  the  components  of  the  FIM  will  be  to  evaluate  a 
representative  term  in  each  block,  followed  by  a  general  expression  for  the  block  itself. 

The  i,j  term  of  the  upper  3x3  ("attitude-attitude")  block  is  found  (A.6)  and  (A.l) 
as 

FIM,,    =    2Re(j>|2d?(Z)Cr1dJ(Z)} 


where  we  have  used  [23] 


(=i 

iAHl 


2ReO    |7>W (OdjtOf  (A.ll) 


2 
and  the  properties  of  Chapter  4 


Re(a)Re(6)  =  -Re(a6*  +  ab)  (A.12) 
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Efw«wT]    =    0 
E[w,wf]    =    Cj6(i-j) 

Using  (A. 11),  the  entire  upper  3x3  block  may  be  written  as 


FIM 


(1-3),  (1-3) 


=    2  Re 


El-nl^^c^KOJ 


2  J>|2D"(/)D(Z) 


i=i 


where 


and 


=    F„ 


D(Z)  =  [d1(0d2(0d3(0] 


6(0=  d1(0d2(/)d3(/) 


(A.13) 
(A-14) 


(A.15) 

(A.16) 
(A.17) 

(A.18) 


(A.19) 


In  the  "attitude-real  gain"  (and  transpose  of  the  "real  gain-attitude" )  block  of  the 
FIM,  the  component  corresponding  to  attitude  term  i  and  gain  j  is 


FIM,,^    =    2Re{7ja//(j)Ci-1di(i)} 

=    2ReUaH(j)di(j)\ 
} 


(A.20) 
(A.21) 


Using  this,  the  entire  row  (column)  corresponding  to  the  ith  attitude  and  the  L  real 
components  of  L  7's  is  found  as 


FIMi,(4,5,...,3+L)    =    2Re{diag{rA"C/-1Ai}} 
=    2Re|diag{rA//AiU 


(A.22) 
(A.23) 


where 


r  = 


7l 

0 

0 

0 

72 

0 

0 

0 

... 

0 

0 

0 

IL  _ 

A  =  [a!  a2  •  •  •  aL] 
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(A.24) 


(A.25) 


A  =  [ai  a2  •  •  •  aL] 
and  the  M  x  L  matrices  Aj  and  A»  are 


A,  =  [^(1)  d,(2)  ■■■di(L)} 


(A.26) 


(A.27) 


Ai  =  [di(l)  4(2)  •  •  •  di(I)J  (A.28) 

Here  diag(x)  forms  a  row  or  column  of  the  diagonal  elements  of  the  square  matrix  x  as 
appropriate.  Mathematically,  this  operation  is  implemented  as 


diag(x)  =  [11  •  •  •  1]  (d  0  I)  (A.29) 

where  I  is  an  identity  matrix.  The  entire  "attitude-real  gain"  (and  transpose  of  the  "real 
gain-attitude")  block  of  the  FIM  is  found  from  the  rows  defined  in  (A.22)  and  (A.23). 


FIM{1_3),(4,5,...,3+L)    =    2  Re 


diagfrA^C^Ai} 

diagjrA^C^Aa} 
diagjrA^Cr'Aa} 
diagjrA^A!} 
diag{rA"A2} 

diag{rA"A3} 

=    Re[Fg,7] 


=    2  Re 


(A.30) 


(A.31) 
(A.32) 
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In  the  "attitude-imaginary  gain"  (and  transpose  of  the  "imaginary  gain-attitude" ) 
block  of  the  FIM,  the  component  corresponding  to  attitude  term  i  and  gain  k  is 

FIMt)L+,+3    =    2Re{^df(A;)Ci-1a,}  (A.33) 

=    2Im{7fca?di(A:)}  (A.34) 

=    2Im{7*i£di(*)}  (A.35) 

using  Re(jx")  —  Im(x). 

Using  the  approach  above,  the  entire  3xL  "attitude-imaginary  gain"  block  is 


FIM 


(1-3),  (4+i,  5+L,  » ,  3+2xL) 


2  1m 


21m 


diagfrA^Cf^j} 

d\ag{rAHCTxA2} 
diagjrA^Cr'Aa} 

diag{rA"Ai} 
diagjrA^Aa} 
diag{rAHA3} 


Im  [F,0 


(A.36) 


(A.37) 


(A.38) 


Using  (A. 14),  the  only  non-zero  terms  of  the  "real  gain-real  gain"  block  of  the  FIM 
are  along  the  diagonal.  The  kth  term  along  the  diagonal  (corresponding  to  the  real  part  of 
7*)  is 


FIM3+*,3+*    =    2Re{afC,-1afc} 
=    2Re{afajt} 


=    2  5fafc 


(A.39) 
(A.40) 
(A.41) 


and  the  entire  "real  gain-real  gain"  block  of  the  FIM  is 


FIM 


1,5,  -,3+i),  (4,5,  -,3+L)      =      A     C,AOl 

=    A"A0l 


£    F, 


(A.42) 

(A.43) 
(A.44) 
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Similarly,  the  "imaginary  gain-imaginary  gain"  block  of  the  FIM  is  also  diagonal.  The 
kth  term  along  the  diagonal  (corresponding  to  the  imaginary  part  of  7*)  is 


FIM 


3+L+k,  3+L+k      — 


2Re{j*jafCi-1afc} 
2Re{afa,} 


-    2afafc 


and  the  entire  "imaginary  gain-imaginary  gain"  block  of  the  FIM  is 


(A.45) 
(A.46) 
(A.47) 


FIM 


(4+L.5+L,  -,3+2xL),  (4+L,  5+L,  -,3+2xL)      =     A     Cj     A©I 

=    AHA0I 
±    F. 


(A.48) 
(A.49) 
(A.50) 


The  remaining  block  is  the  "real  gain-imaginary  gain"  block.  From  (A.4)  and  (A. 5), 
the  FIM  component  corresponding  to  real  gain  i  and  imaginary  gain  A;  is 


FIM 


3+i,  3+L+k 


E  [4  Re  {af  GTh*}  Im  {af  Cflwfc}] 

E  [af  CT1w<wJ,CJ-r«J  +  (•? cr^cr'-DI 


0V«,fc  =  1,2, 


(A.51) 
(A.52) 
(A.53) 


Using  (A.17),  (A.32),  (A.38),  (A.44),  (A.53),  and  (A.50),  the  FIM  may  be  written  as 


FIM  = 


F 

r9 

Re[F,,7] 

Im  [F,, 

Re  [F,,7]T 

F7 

0 

.  Im  [F9,7]T 

0 

F 

(A.54) 


Of  interest  in  this  problem  is  the  CRB  for  the  attitude;  the  complex  gain  terms  are 
essentially  nuisance  parameters.  The  CRB  for  the  desired  attitude  parameters,  PCT(9),  is 
found  from  the  upper  3x3  block  of  the  inverse  of  the  FIM.  Using  the  inverse  identity  of 
Liet  al.  [23], 

PCT(9)-1  -  Re  {Fg  -  F^F^F^}  (A.55) 
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Now  examine  the  terms  of  P"1^).  Since  F7  is  a  diagonal  matrix,  (F7)_1  is  also  diagonal. 
Using  this,  the  i,  j  element  of  F^F^F^,  is  found  to  be 

a.H(l)C^dj(l)  (A.56) 

L 

=    2^|7,|2df(/)P5(()dJ(/)  (A.58) 


(=1 


The  entire  3x3  matrix  of  F^F^F*  is 


F^F^Fj,  =  2  <T  |7(|2Df  (J)P«<,)D(Z)  (A.59) 
Finally,  recall  that  F9  =  2  Re  fcL=1  \ji\2DH {l)B(l)V  so  that 

P"1^)    -    2  Re  j  ^|7;|2{DH(/)D(/)-DH(/)PaWD(/)}|  (A.60) 

=    2Re(x:|7i|2D//(/)Pi(()D(Z)}  (A.61) 


(=i 


APPENDIX  B 
THE  INTERFERENCE  COVARIANCE  MATRIX 

B.l     Introduction 

One  of  the  fundamental  requirements  in  the  development  of  the  maximum  likelihood 
methods  of  this  work  is  the  knowledge  of  the  statistics  of  the  underlying  interference. 
In  this  chapter  two  important  topics  are  presented,  the  theoretical  calculation  of  the 
interference  and  the  estimation  of  the  interference  for  a  practical  system  .  A  two  sensor 
system  and  a  simplified  receiver  chain  are  used  for  the  analytical  calculation,  as  they 
provides  a  clear  illustration  of  the  underlying  phenomena.  This  two  sensor  system  is  used 
to  calculate  the  space-satellite  covariance  matrix  defined  in  Chapter  4.  Once  the  model  for 
the  interference  has  been  developed,  methods  of  estimating  the  interference  are  discussed. 
In  particular,  a  method  specific  to  this  application  is  developed,  which  may  prove  useful 
from  a  computational  standpoint. 

In  this  chapter  we  remove  the  constraint  of  Chapter  4  that  the  interference  is  tem- 
porally uncorr elated  and  instead  calculate  the  various  interference  statistics  for  a  jammer 
with  a  band-limited  power  spectral  density.  The  band-limited  spectra,  the  derivations  of 
the  interference  statistics,  and  the  methods  of  estimating  the  interference  developed  in  this 
chapter  are  implemented  in  the  simulation  used  to  produce  the  results  of  Chapter  8. 

In  the  derivations  that  follow,  accurate  representations  of  the  spatial  covariance 
matrices  are  developed  from  the  jammer  temporal  correlation  functions.  This  is  to  ensure 
accuracy  in  the  equations  as  the  spacing  between  sensors  increases.  As  the  array  spacing 
increases,  the  jammer  signals  experience  decorrelation  between  the  various  sensors.  For 
small  spacing,  the  sensor  to  sensor  correlation  is  adequately  approximated  by  a  phase 
difference  between  sensors  (i.e.  no  amplitude  decorrelation).  In  this  case  an  adequate 
model  for  the  jammer  is  simply  the  array  response  vector,  which  is  composed  of  these 
phase  differences.  However,  as  the  distance  between  sensors  increases  or  the  jammer 
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bandwidth  increases,  the  amplitude  of  the  sensor  to  sensor  cross  correlation  is  attenuated 
in  addition  to  the  phase  shift.  This  amplitude  decrease  in  fact  degrades  performance  over 
the  phase-only  case. 1    Since  Chapter  7  focuses  on  wider  baseline  lengths  than  typical 
(for  array  processing)  half  wavelength  spacing,  accurately  representing  the  performance 
of  the  attitude  determination  algorithms  of  this  work  required  a  more  detailed  jammer 
decorr elation  approach  than  the  phase-only  method. 

Nomenclature 

This  chapter  uses  a  considerable  amount  of  nomenclature.  To  remove  confusion,  the 
important  terms  and  their  definitions  are  listed  below. 

At  the  sampler,  the  interference  terms  from  Chapter  4  are  composed  of  contributions 
from  thermal  noise  and  jammers. 

n[A;]    =    kth.  sample  of  interference 

nt[k]    =    Thermal  Noise  Contribution  to  A;th  sample  of  interference 
iij[A;]    =    Jammer  Contribution  to  kth  sample  of  interference 

The  spatial  interference  covariance  is  formed  from  samples  taken  at  the  A/D,  i.e. 
before  despreading.  This  is  composed  of  contributions  from  the  jammer  and  thermal  noise. 

Rs    =    Spatial  Covariance 

Rs,    =    Jammer  Contribution  to  Spatial  Covariance 
Rsf    =    Thermal  Noise  Contribution  to  Spatial  Covariance 

(B.l) 


1  In  the  limiting  case  of  no  sensor  to  sensor  correlation,  the  jammer  power  appears  spa- 
tially white,  like  thermal  noise.  Here  performance  may  be  significantly  reduced,  as  the 
situation  may  be  viewed  as  an  unjammed  case  with  the  thermal  noise  power  equal  to  the 
sum  of  the  actual  thermal  noise  power  and  the  (possibly  very  large)  jammer  power. 
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After  demodulation,  the  interference  is  composed  of  contributions  from  the  jammer, 
thermal  noise,  and  multiple  access  interference  (MAI)  from  other  satellites. 


w, 

Wlt 

Wlm 


M  x  1  Interference  Vector  in  the  /th  Satellite  Channel 
Jammer  Contribution  to  the  Ith  Interference  Vector 
Thermal  Noise  Contribution  to  the  Ith  Interference  Vector 
MAI  Contribution  to  the  Ith  Interference  Vector 


The  spatial  -satellite  interference  covariance  is  formed  from  samples  taken  after 
despreading.  This  is  composed  of  contributions  from  the  jammer,  thermal  noise,  and 
multiple  access  interference  (MAI)  from  other  satellites. 

Rss  =  Spatial  -  Satellite  Covariance 

Rssj  sa  Jammer  Contribution  to  Spatial  -  Satellite  Covariance 

Rsst  =  Thermal  Noise  Contribution  to  Spatial  -  Satellite  Covariance 

RsSm  =  MAI  Contribution  to  Spatial  -  Satellite  Covariance 

Rss[/m]  =  The  l,m  M  x  M  Block  of  Spatial  -  Satellite  Covariance 


(B.2) 
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Figure  B.l:  The  receiver  model  used  in  this  dissertation.  The  output  of  the  chip  matched 
filter  is  oversampled  by  a  factor  of  JVC)  (Nc  is  4  in  this  figure)  and  these  samples  are  deci- 
mated into  sequences  composed  of  1  sample  per  chip.  The  decimated  sequences  are  then 
despread  using  the  sequency  A(i)    i    =    1,2, ..  .N  to  form  early,  punctual,  late,  and  other 
temporal  gates. 
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The  general  receiver  design  of  Figure  B.2  is  similar  to  the  DS-SS  receiver  described  in 
[45].  This  system  employs  an  anti-aliasing  filter  matched  to  the  chip  waveform,  and  then 
oversamples  the  output  of  this  filter.  The  samples  are  then  decimated  and  each  decimated 
sequence  is  individually  de-spread  to  provide  a  punctual  output  and  possibly  several  other 
outputs  delayed  and  advanced  in  time  as  required  for  acquisition  and  tracking.  With  this 
method,  only  one  sample  per  chip  is  used  in  any  correlation  (despreading)  function,  and 
up  to  Nc  correlators  could  be  implemented  in  the  time  span  of  one  chip. 

The  spatial  interference  covariance  matrix  is  composed  of  contributions  from  all  jam- 
mers and  thermal  noise.  A  few  reasonable  assumptions  allow  straightforward  calculation 
of  this  covariance.  First,  we  assume  that  all  jamming  sources  are  zero  mean,  circularly 
symmetric,  wide-sense  stationary  complex  Gaussian  random  processes,  and  that  each 
jammer's  waveform  is  independent  from  those  of  all  other  jammers.  Second,  we  assume 
that  the  physical  nature  of  the  thermal  noise  causes  it  to  be  uncorrelated  from  channel 
to  channel  (spatially  white),  and  that  its  power  spectral  density  is  uniform  over  the  fre- 
quencies of  interest  (temporally  white).  Finally,  we  assume  that  all  jammer  waveforms  are 
uncorrelated  from  all  satellite  waveforms. 

B.2     Derivation  of  Interference  Statistics 
Space-Satellite  Interference 

The  required  statistics  of  all  interference  sources  are  captured  in  the  space-satellite 
covariance  matrix  RsS.  To  review  from  Chapter  4,  this  is  composed  of  the  expected 
value  of  the  outer  product  of  the  space-satellite  interference  vector,  W,  a  ML  x  1  vector 
containing  the  interference  in  each  satellite  channel. 

Rss  =  E  [WWH]  (B.3) 
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r              -i 

Wi 

w  = 

w2 

wL 

(B.4) 


For  this  work,  the  interference  is  comprised  of  contributions  from  three  types  of  source: 
jammers,  interfering  satellites  (i.e.  multiple  access  interference),  and  thermal  noise. 


W|  =  wij  +  wIt  +  w!r 


(B.5) 


This  can  be  seen  from  equation  (4.37),  the  interference  in  each  satellite  channel  is  com- 
posed of  the  thermal  noise  and  jammers,  and  the  multiple  access  interference  (MAI). 

L 


1 

Wj  =  - 

ei 


[n[ti],n[*2],...,n[*tf]]+    ^    lPap(q)yf 

p=i,  p#( 


Yt 


H 


(B.6) 


since  the  n[k]  terms  were  defined  in  Chapter  4  to  contain  contributions  from  thermal  noise 
and  jammers,  i.e. 

n[k]  =  nj[k]  +  nt[k]  (B.7) 

As  is  the  assumption  throughout  this  work,  each  of  these  are  independent  from  one 
another,  and  therefore  the  total  interference  realization  and  its  covariance  are  simply  the 
sum  of  each  of  the  individual  realizations  and  covariances,  respectively. 


—  -H-ssi  T  -Ttsst    i    -t*-ssr 


(B.8) 


Since  the  jammer  and  thermal  noise  components  of  the  space-satellite  covariance  matrix 
are  found  from  the  spatial  covariance  Rs,  we  begin  by  examining  the  components  of 
spatial  covariance. 

Jammer  Contribution  to  R. 


This  section  calculates  the  interference  covariance  for  one  jammer.  Since  each  jammer 
(if  more  than  one  exist)  is  independent  from  all  others,  their  contributions  would  all  be 
calculated  independently  using  the  method  described  in  this  section,  and  simply  summed. 
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The  receiver  model  described  earlier  in  this  chapter  consists  of  a  mixer,  which 
downconverts  the  random  process  to  baseband,  a  matched  filter,  and  a  sampler  (A/D  con- 
verter). Following  the  sampler,  the  direct  sequence  spread  spectrum  (DS-SS)  modulation 
is  removed,  and  the  de-modulated  samples  are  summed,  resulting  in  an  output  data  vector 
wjj,  I  =  1,. . .  ,L  which  is  used  for  the  attitude  determination.  Repeating  this  process  with 
various  time-shifted  replicas  of  the  spreading  sequence  generates  the  appropriate  data  for 
code  tracking  (e.g.  late  gate,  early  gate,  etc.). 

We  define  the  random  process  corresponding  to  the  jammer  signal  that  we  will 
evaluate  as  x(t).  The  shape  of  power  spectral  density  (PSD)  is  assumed  to  be  flat  across 
the  frequency  band  fc  —  B  <  f  <  fc  +  B,  where  fc  is  the  GPS  center  (carrier)  frequency. 
The  autocorrelation  function  corresponding  to  this  flat  spectrum  has  the  familiar  sine 
shape  in  amplitude,  and  a  phase  related  to  the  carrier  frequency  and  time  delay.  Px  is  the 
power  (variance)  of  the  process. 


'c  -B        f0     fc  +B 


frequency 


Figure  B.2:  One  sided  jammer  PSD 


Rxx(t)    =    E[x*(t)x(t  +  T)} 

=    p  cj»cts™(^Bt) 

tt2Bt 
=    Pxej^Tsinc(2BT) 


(B.9) 
(B.10) 


where 


sinc(a)  = 


sin(7ra) 


na 


(B.ll) 


The  two  sensors  receive  this  random  process,  with  the  second  sensor  receiving  the 
process  delayed  by  an  amount  of  time  A,  which  is  related  to  the  direction  of  arrival  of  the 
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source  and  the  physical  arrangement  of  the  two  sensors.  Denote  the  two  signals  received  in 
the  two  sensors  as  s\(t)  and  S2(t). 

8t{t)    =    x(t)  (B.12) 

s2(t)    =    x(t-A)  (B.13) 

The  receiver  will  downconvert  these  bandpass  processes  to  baseband  random  pro- 
cesses. We  denote  the  baseband  processes  in  each  channel  corresponding  to  the  downcon- 
verted  waveform  x(t)  as  Si(t)  and  52 (£). 

Si(t)    =  «i(e)e-*^  (B.14) 

=  xtyy-P**  (B.15) 

S2(t)    =  s2(i)e-jWct  (B.16) 

=  x{t-A)e-j^t  (B.17) 

The  process  Si  (t)  is  a  zero  mean  normal  wide  sense  stationary  random  process  with 
autocorrelation  function  RsltSi(T~)- 

=  E[x*  (t)e+JWctx(t  +  r)e-^(t+T)] 

=  E[x*(t)x{t  +  T)}e-ju,^T) 

=  PJSetu'Tmac(2BT)e-*^r) 

=  PlSinc(25r)  (B.18) 

Since  x(t)  is  wide  sense  stationary,  then  52(^)  is  wide  sense  stationary  as  well. 
In  addition,  Rsia(t)  and  Rs2,s2(T)  are  identical,  which  can  be  easily  seen  by  letting 
y(t)  =  x(t  -  A)  and  observing  that  the  statistics  of  y  and  x  are  identical. 

RSi,Si(t)  =  RS2,S2(t)  (B.19) 
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The  cross  correlation  function  is  found  in  a  similar  manner. 


rftA(T)    =  E[Sl(t)S2(t  +  r)} 

=  E[z*(*)e+M;fx(t  +  r  -  A)e-^c(t+T)] 

=  E[i*(*)x(i  +  r-A)]e-^T 

=  Rxx(t  -  A)e-j"°T 

=  Pxsinc  (2B(r  -  A))  e**M)e-WT) 

=  Pxsinc(25(r-A))e-^A  (B.20) 

Two  observations  may  be  made  in  comparing  the  cross  correlation  term  Rsus2  to  the 
autocorrelation  term  (either  i?Si,Si  or  Rs2,s2)-  First,  the  amplitude  of  the  cross  correlation 
term  at  r  =  0  is  reduced  by  an  amount  corresponding  to  the  time  delay  between  reception 
at  the  sensors  and  the  process  bandwidth.  Second,  a  carrier  frequency  phase  shift  is 
imposed. 

Following  the  downconverter,  the  signal  in  each  channel  is  matched  filtered,  sampled, 
and  decimated  by  a  factor  of  JVC,  producing  two  discrete  time  sequences,  one  per  channel. 
The  matched  filter  for  the  BPSK  type  GPS  signals  is  simply  an  ideal  integrator.  In 
practice,  however,  a  small  SNR  penalty  could  be  paid  for  deviating  from  the  matched 
filter  to  one  with  a  sharper  transition  region,  greater  stop  band  rejection,  or  easier 
implementation  (see,  for  example,  Chapter  10  of  [46]).  We  denote  these  sequences  as  Si[k] 
and  S2[k],  where  k  is  the  time  (sample)  index.  The  integration  time  the  chip  period  Tc. 

rkTc 


PK1C 

Si[k]  =  /  S!(t)dt  (B.21) 

J(k-1)TC 
rkTc 

S2[k]  =  /  S2(t)dt  (B.22) 

J(k-l)Tc 


Notice  that  S\[k\  and  S2[k]  are  zero  mean  normal  wide  sense  stationary  random  sequences, 
as  shown  below  for  by  its  mean, 

rkTc 

E[5i[*}]    =    /         E[Sx(t)]dt 

J(k-1)TC 

=    0  (B.23) 


and  its  autocorrelation: 


rslA[J]    =    E[SJ[fc]Si[*  +  Z]] 

rkTc  r(k+l)Tc 


=  /  E[5J(r)5i(A)]rfTdA 

J(k-1)TC  J(k+l-l)Tc 
rkTc  r(k+l)Tc 

=  /  /2SliSl(A-T)(iTdA 

J(k-1)TC  J{k+l-l)Tc 
rTc     Ml+l)Tc 

=  Pxsmc(2B(X  -  T))drdX 

J0      JlTc 


149 


(B.24) 


The  cross  correlation  function  between  the  discrete  sequences  is  similarly  found  by  direct 
evaluation. 


r«A[(l         E[Sl[k}S2[k  +  l}} 

rkTc  r(k+l)Tc 

E[Sl(T)Si(\))<ird\ 

'(lt-l)Tc  J(k+l-l)Tc 
rkTc  r(k+l)Tc 

Rsi&tft  -  T)drdX 

l(k-\)Tc  J{k+l-\)Tc 
rTc     r(l+l)Tc 

=  Pxsmc(2B(\-T-A))e-j"°AdTd\ 

J0      JlTc 


fkTc  A 

J(k-iyr,  J{k 

rkTc  r( 

Jtk-l)Tc  J(k 


(B.25) 


The  sequences  5X  [A;]  and  S2  [k]  form  the  vector  of  jammer  contributions  to  interference 

at  the  A/D,  n.j[k}. 

""  S1[k] 


Uj[k]  = 


S»[k] 


(B.26) 


The  jammer  covariance  at  the  A/D,  KSj  is  obtained  from  Uj[k]  as 


RSJ  E  [n,[A;]nf  [k]] 

rSuSl[0]    rM|Dj 
r*Sl,sM    r52,S2[0] 

This  can  be  written  concisely  to  show  the  two  cross  correlation  properties  of  phase 
change  and  amplitude  reduction. 


(B.27) 
(B.28) 


R, 


*: x*< rr 


1  (1  -  a)  e-l°°* 

(I -a)  ej"'A  1 


(B.29) 
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where 


rkTc          rkTc 

Panto    =                  /           FIsinc(2B(A- 

J{k-\)TC  J(k-1)TC 

-  r))drdX 

(B.30) 

rkTc           rkTc 

Pcross    =                  /           PIsinc(25(A- 

J(k-\)TC  J(k-1)TC 

-  T  -  A))dTd\ 

(B.31) 

&      —      1         *  cross/ •» auto 

(B.32) 

Approximation  for  Narrow-Band  Jammers 

Consider  the  diagonal  terms  of  (B.29)  for  narrow-band  processes.  If  the  time  band- 
width product  BTC  is  small,  then  the  sinc(£?Tc)  function  is  essentially  constant  and  equal 
to  1.  In  this  case,  nearly  all  of  the  power  in  x(t)  passes  through  the  matched  filter. 

rkTc  rkTc 

rsusM    =  /  Pxsmc{2B(X-T))dTd\ 

J(k-1)TC  J(k-1)TC 
rTc     rTc 

Prdrd\ 


Jo     Jo 


w    PxTl  (B.33) 

A  similar  narrow-band  analysis  may  be  performed  for  the  off-diagonal  terms  of  (B.29)  by 
using  the  cross-correlation  term.  However,  in  this  case  the  additional  delay  corresponding 
due  to  the  difference  in  signal  arrival  times  between  the  two  sensors  must  be  considered. 
If  the  time  bandwidth  product  B(TC  +  A)  is  small,  then  the  sine  function  is  essentially 
constant  and  the  cross  correlation  value  is  equal  to  the  autocorrelation  value  multiplied  by 
a  complex  exponential  phase  shift. 

rkTc  rkTc 

r5,,sa(0)    =     /  /  Pxsinc(2B(X-T-A))e-juJ°AdTdX 

J(k-\)TC  J{k-\)TC 
Tc     fTe 

Pxe-j^AdTdX 


II 

Jo     Jo 


ps    PxT^e~iUeA 


=    e-JulcArSuSl(0)  (B.34) 
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For  the  small  time  bandwidth  cases  described  above,  Pauto  and  Pcross  are  essentially  equal 
to  Tc2,  and  Rs  [0]  is  now  the  rank  1  matrix. 


Rs .  ea 


PXT'C 


1        e-**A 


(B.35) 


The  eigenvector  of  this  matrix  corresponding  to  the  non-zero  eigenvalue  is  the  array 
response  vector  (i.e.  spatial  steering  vector)  associated  with  the  direction  of  arrival  of  this 
interference  source.  In  many  array  processing  applications  this  steering  vector  model  is 
adequate.  However,  since  this  dissertation  explores  large  array  spacing,  it  uses  the  more 
accurate  amplitude  decorrelation  version. 

Thermal  Noise  Contribution  to  R, 


The  non-zero  temperature  of  physical  devices  causes  motion  of  electrons  (i.e.  current) 
inside  the  material,  resulting  in  a  random  voltage  at  the  terminals  of  the  device  [47].  This 
thermal  noise  is  an  important  factor  in  system  performance.  The  common  "Additive 
White  Gaussian  Noise"  (AWGN)  model  for  thermal  noise  considers  the  noise  to  be 
spectrally  white,  normally  distributed,  and  independent  of  all  other  signals  in  the  receiver. 
Because  of  the  physical  differences  between  receiver  channels,  the  noise  one  receiver 
channel  is  assumed  independent  from  that  in  all  the  other  channels. 

Define  ntj  (t)  as  the  thermal  noise  waveform  appearing  in  the  2th  hardware  channel. 
After  matched  filtering  and  sampling,  the  resulting  random  sequence  ntJfc],  i  —  1,2  forms 
the  thermal  noise  contribution  to  the  interference  vector  nt[k}. 


fkTc 

"*[*]=  / 

J(k-l) 


(fc-l)rc 


nt[k]  = 


nti(t)dt        i  =  l,2 

nt2[k] 


(B.36) 


(B.37) 
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Using  the  assumptions  imposed  on  nt[k]  the  thermal  noise  covariance  at  the  A/D,  Rst 
is  the  diagonal  matrix  Rgt 


(B.38) 
(B.39) 


Rsf  E  [nt[k}n?lk}] 

Ptjc         0 

0       Pt2Tc 

where  Ptl  and  Pt2  are  the  power  spectral  density  level  of  the  thermal  noise  in  channels  1 
and  2,  respectively. 

Spreading  Sequence  Model 

The  demodulation  waveform  yi  is  the  complex  sequence  composed  of  the  spreading 
sequence  and  estimated  Doppler  to  the  /th  satellite.  Although  this  sequence  is  determin- 
istic (it  is,  of  course,  known  to  the  receiver)  it  is  useful  to  model  it  as  a  random  process 
when  the  spreading  gain  is  large  (i.e.  the  number  of  chips  integrated  is  large).  The  model 
for  the  random  process  uses  the  common  approach  of  considering  each  of  the  terms  in  the 
sequence  to  be  independent,  identically  distributed  variables,  giving  rise  to  the  following 
properties: 

E[yi[n}y;[n-k}}  =  6{k]  (B.40) 

Etei[n]Vt[k]}  =  Ql?i,  Vn,*  (B.41) 

where 

1    k  =  0 

0    k^O 

These  above  relationships  flow  from  the  asymptotic  properties  of  the  GPS  P(Y)  code 

spreading  sequences: 


S[k] 


(B.42) 


1    N 
iim  7  2J  KMtfl"  -  k]  =  °         k  ^  ° 


71=1 


A' 


,fec  7  S  y'Ny?  [n  -  *]  «  0         V  A: 

'  n=l 


(B.43) 


(B.44) 
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Where  e,  the  energy  in  the  z'th  sequence,  is  defined  as 

N 


=  Z>HtfM  (B-45) 


n=l 

Jammer  Contribution  to  RgS 

To  calculate  the  jammer  contribution  to  the  space-satellite  interference  covariance 
matrix,  we  partition  RsSj  into  M  x  M  blocks,  and  consider  blocks  on  the  diagonal 
and  blocks  off  the  diagonal  separately.  First  consider  the  ith  diagonal  block,  which  is 
calculated  from  w^,  the  jammer  waveforms  despread  using  the  zth  despreading  sequence: 

1    - 
«fc-7X>M*fW  (B.46) 

Using  this,  the  ith  diagonal  block  of  the  jammer  contribution  to  Rss  is  found  from 


YUsj[ii]    =    E[wJiWjf]  (B.47) 

"    3'EEBhW^fWIi*4*WI  (B-48) 


1    a=l  0=1 

Using  the  statistical  properties  of  the  spreading  sequence  defined  in  equation  (B.40)  (i.e. 


1      N      N 

*W«1    *    ?£EE[*^]E|»;[aW]  (B.49) 


treating  jft  as  sequence  of  independent  random  variables) ,  this  evaluates  to 

N      N 

1   a=l  0=1 

N 

=  sER-ibiW!2  (B.50) 

*    a=l 

=    ±RS]  (B.51) 

■ 

This  is  an  important  result:  the  diagonal  blocks  of  the  space-satellite  covariance  matrix 

are  all  equal  to  the  jammer  spatial  covariance  matrix. 

In  a  similar  fashion,  the  off  diagonal  blocks  may  be  calculated  as  well.  Consider  the 

i  —  j  block 

1  1    N    N 
*~M  =  77  E  E  E  [»iWnf  Wl  E  WW*WI  (B-52) 

Q=l  /3=1 
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Again  using  the  properties  of  the  spreading  sequence  (specifically  equation  (B.41)),  the 
off-diagonal  terms  evaluate  to  zero. 

Rss^]  =  0  (B.53) 

Combining  the  diagonal  and  off-diagonal  terms  produces  the  total  contribution  of  the 
jammer  to  the  space  satellite  interference  covariance  matrix.  Since  the  energies  in  each 
of  the  spreading  sequences  is  the  same,  we  may  drop  the  satellite  subscript  on  the  energy, 
and  write  the  jammer  contribution  as 


Rs 


sj 


R*,     o 


o 


0 

R-sj 

0 

0 

0 

... 

0 

0 

0 

R-sj 

(B.54) 


Thermal  Noise  Contribution  to  Rss 

Using  the  same  methodology  as  for  the  jammer,  the  contribution  from  the  thermal 

noise  is  found  to  be 

Rst      0      ■■•      0 

0      Rs(     0       0 

o     •••    ••.     0 

0       •••      0     YUt 

and  since  Rgt  is  a  diagonal  matrix,  RsSt  is  as  well. 


R-ss 


(B.55) 


MAI  Contribution  to  RsS 

The  third  type  of  interference  considered  is  that  from  satellites  other  than  the  desired 
satellite  in  each  satellite  channel,  i.e.  the  multiple  access  interference  (MAI).  Conceptually, 
if  there  are  L  satellites  being  tracked  and  used  for  attitude  estimation,  then  in  the  first 
satellite  channel  satellites  2  -  L  are  interferers.  Similarly,  in  the  second  channel,  satellites 
1  and  3  -  L  axe  interferers,  and  so  on.  The  MAI  in  the  Ith.  satellite  channel,  wlm,  is 
composed  of  contributions  from  all  satellites  but  the  Zth: 

L 


wiro  =  - 
Si 


J2  v»p(9)y* 

p=l,pjil 


yf 


(B.56) 
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To  calculate  the  contribution  of  MAI  to  the  space-satellite  interference  covariance 
matrix,  we  partition  RgS  into  M  x  M  blocks,  and  consider  blocks  on  the  diagonal  and 
blocks  off  the  diagonal  separately. 


Wm  = 


WLm 


Rssm  =  E  [WmW£] 
The  zth  diagonal  block  of  RSSm  is  found  to  be: 

,[H]    =    E[w,TOw,£] 


1 


=    -^E 


Y     Y   7»7,*ai(5)af (ftsWiW? 


'      \_i=\  i^i  j=n^i 


1 


=    ±   £    l7*l^(5)af  (5) 
'  »=i  <#( 

The  /  —  p  off-diagonal  block  is  calculated  to  be: 


M  = 


i 


E 


SlEp 
=     0 


(B.57) 


(B.58) 

(B.59) 
(B.60) 

(B.61) 


PmJ 

(B.62) 

'      L             L 

Y     Y   7i7J*ai(g)af(g)yiyi*ypy* 

(B.63) 

i=l  i^i  j=l  i^p 

(B.64) 

Combining  these  diagonal  and  (zero)  off-diagonal  terms,  the  total  MAI  contribution 
to  the  space-satellite  interference  covariance  matrix  is  a  block  diagonal  matrix: 


"«sm  — 


Rssm  [11]                0 

0 

0                Rssm  [22] 

0 

0 

0 

0 

0 

0 

Rssm  [LL\ 

(B.65) 
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An  important  observation  is  that  since  the  satellites  that  compose  the  MAI  in  each 
satellite  channel  are  different,  each  of  the  blocks  that  make  up  Rssm  are  different  as  well. 
It  is  for  this  reason  that  in  the  derivations  of  Chapter  4,  the  covariances  are  indexed  by 
the  satellite  channel  subscript. 

B.3     Estimation  of  Interference  Statistics 
It  is  very  important  in  estimation  of  the  statistics  of  the  interference  to  ensure  that 
the  training  data  for  the  estimation  are  free  of  the  desired  signal.  If  the  desired  signals 
exist  in  the  training  data,  the  estimate  of  the  desired  parameter  (attitude  or  direction  in 
this  work)  will  be  degraded.  For  a  pulsed  radar,  estimation  of  the  interference  may  be 
performed  during  the  portion  of  the  pulse  repetition  interval  corresponding  to  where  the 
desired  target  is  not  expected  to  be  (i.e.  ranges  away  from  the  target).  However,  for  a 
communication  system  where  the  communication  waveforms  are  present  at  all  times,  other 
methods  must  be  employed.  In  this  section  we  briefly  discuss  methods  to  estimate  the 
space-satellite  interference  covariance  matrix  Has. 

Estimation  of  R^Jm] 


Knowledge  of  the  time  waveforms  of  the  signals  provides  an  advantage  in  estimating 
the  interference  covariance.  Since  for  any  satellite  channel  the  desired  (matched)  satellite 
waveform  lies  in  1  dimension  of  an  N  dimensional  space  containing  all  possible  length- N 
BPSK  signals,  the  desired  signal  may  be  removed  for  interference  estimation  by  projecting 
the  received  data  onto  a  subspace  of  this  N  dimensional  space  that  is  orthogonal  to  the 
satellite  time  signals.  Examples  of  this  method  of  time-series  projection  may  be  found  in 
Vaidyanathan  et  al.  [48]  for  CDMA  communication  system  with  an  antenna  array  that 
does  not  make  use  of  the  LOS  to  the  source,  in  Jones  and  Wikert  [49]  for  a  robust  system 
that  incorporates  both  LOS  information  (in  a  generalized  side-lobe  canceller  architecture) 
and  knowledge  of  the  desired  waveform,  and  in  Li  et  al.  [22, 23]  for  a  direction  finding 
application. 

Recall  that  the  M  x  N  data  matrix  X  of  equation  (4.25)  contains  the  data  samples 
from  each  sensor  for  N  time  snapshots.  Using  this,  an  estimate  of  the  tth  diagonal  block 
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of  the  space-satellite  covariance  matrix  can  be  found  by  "projecting  out"  (i.e.  projecting 
onto  a  space  orthogonal  to)  the  ith  satellite  waveform: 

RsJm] 


1 

\xXH         1//(Xyf)(yiXw)l 

1 

XINXH  -^XPyiX" 

=  -1 

Ns 

X(I^-Pi)Xff] 

TK  [XPv-X"] 


Ne 


(B.66) 


where 


P*    = 


Yi   Yi 


YiY 


B 


Py.     =    ^"J^P* 


(B.67) 
(B.68) 


i=i 


GPS  Receiver  Considerations 


Recall  from  Section  2.2  that  the  signal  level  of  the  satellite  waveforms  at  the  sampler 
is  much  less  than  that  of  even  the  thermal  noise.  This  requires  that  many  chips  be  de- 
spread  and  integrated  to  provide  appreciable  SNR,  and  as  discussed  earlier  this  function 
is  typically  implemented  in  hardware.  Therefore,  to  calculate  covariance  snapshots  using 
the  entire  data  matrix,  as  described  in  (B.66),  would  likely  be  computationally  prohibitive 
and  unnecessary.  In  addition,  it  would  be  desirable  to  explore  the  potential  to  simplify 
calculations  by  using  a  single  M  x  M  interference  covariance  matrix,  instead  of  L  matrices. 

To  explore  these  two  desires,  consider  augmenting  the  receiver  by  an  additional 
satellite  (vector)  demodulation  channel,  and  denoting  this  as  channel  0.  As  in  (4.35),  the 
output  of  this  channel,  Uo,  is  composed  of  the  sum  of  satellite  and  interference  projections 
onto  the  demodulation  waveform,  y0. 


U0  =  Z0  +  W0 


(B.69) 


Now  if  the  receiver  is  designed  to  collect  the  output  of  channel  0  at  a  rate  much  faster 
than  the  other  channels  (i.e.  much  faster  than  the  attitude  update  rate),  yet  much  slower 


158 

than  the  system  sample  rate2  ,  and  if  we  choose  the  channel  0  waveform  properly,  then  we 
may  estimate  the  interference  covariance  matrix  from  the  samples  of  Uo-  This  will  produce 
the  same  covariance  matrix  for  each  satellite  channel.  If  JVj  is  the  number  of  snapshots  of 
channel  0  available,  the  the  covariance  estimate  is  found  as 

R«[«]  =  J^uo[n]        t  =  l,2,...I  (B.70) 

n=l 

In  order  to  employ  this  method,  an  appropriate  demodulation  waveform  for  channel  0 
must  be  chosen.  Assume  that  a  waveform  yx  exists  that  satisfies  the  two  properties  of  the 
spreading  sequence  model  described  in  equations  (B.40)  and  (B.41).  That  is, 

E[y1[n}y*1[n-k}}=5[k]  (B.71) 

and 

E[y±[n]y;[k]]=OljLi,Vn,k  (B.72) 

then  by  using  yx  is  the  demodulation  sequence,  K^ [it]  is  an  approximation  to  each  of 
the  diagonal  blocks  of  R^n].  The  difference  is  that  since  all  satellites  are  present  in  this 
data,  the  desired  satellite  is  present  in  the  estimated  interference  covariance  matrix  in  each 
channel,  where  in  the  true  covariance  matrix  the  desired  signal  is  absent.  However,  the 
desired  satellite's  contribution  is  attenuated  significantly,  since  its  projection  onto  y±  is 
small.  In  fact,  its  amplitude  is  attenuated  to  the  level  of  the  MAI,  which  is  very  small  for 
the  GPS  P(Y)  waveform  in  general. 

To  summarize,  employing  the  above  covariance  estimation  approach,  two  important 
gains  are  made  in  computational  efficiency.  First,  the  number  of  outer  products  is  reduced 
by  using  multiple  snapshots  of  despread  data  to  estimate  the  covariance,  as  opposed  to 


2  For  example,  if  the  attitude  update  rate  is  100Hz,  and  the  GPS  waveforms  are  sam- 
pled at  10MHz,  then  the  output  of  channel  0  could  be  at  10kHz.  This  would  provide  100 
samples  from  which  to  estimate  covariance,  and  the  spreading  gain  would  be  1000,  suffi- 
cient to  satisfy  the  asymptotic  assumptions  made  earlier.  In  general,  the  number  of  snap- 
shots from  which  to  estimate  covariance  must  be  several  times  greater  than  the  number  of 
sensors. 
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the  spread  data  at  the  sampler.  Second,  by  paying  a  small  price  in  accuracy,  a  single 
covariance  matrix  may  be  used  for  all  satellite  channels,  reducing  by  a  factor  of  L  the 
computational  requirements  for  estimating  the  covariance. 


APPENDIX  C 
RESOLUTION  OF  THE  ATTITUDE  AMBIGUITY  BY  TWO  SATELLITES 

In  this  appendix  we  review  the  attitude  ambiguity  and  show  that  only  two  sources 
are  required  to  eliminate  the  ambiguity  in  the  array  response  vector.  We  begin  by  defining 
the  terms  used  in  this  chapter.  Let  vx  and  v2  be  the  LOS  vectors  to  two  sources  such 
that  vx  7^  vi.  Also  let  qi(wi)  represent  the  locus  of  possible  attitudes  that  produce  the 
same  array  response  vector  from  source  vi,  and  let  c^i^)  represent  the  locus  of  possible 
attitudes  that  produce  the  same  array  response  vector  from  source  V2.  Finally,  define  the 
true  attitude  as  q,  and  partition  its  components  into  the  3x1  vector  z  and  the  scalar  s. 


(C.l) 


As  presented  in  Chapter  4,  the  locus  of  possible  attitudes  from  either  of  the  two 
sources  may  be  written  as  the  quaternion  multiplication 


qi(wi)  =  qi(u>i)*q        *  =  1,2 


(C.2) 


where 


q»(^i) 


t  =  l,2 


(C.3) 


sin(if)vi 
cos(f) 

Using  these  definitions,  we  present  the  following  theorem. 
Theorem  4   The  attitude  ambiguity  in  the  array  response  vector  for  one  satellite  intersects 
the  ambiguity  in  the  array  response  vector  for  a  satellite  at  a  different  direction  in  a  single 
attitude,  the  true  attitude  of  the  antenna  array.   That  is, 


qi(a>i)nq2(cj2)  =q 


(C.4) 
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Proof:     Begin  by  noting  that  when  the  two  attitude  quaternions,  qi(wi)  and  q2(tt>2), 
are  equal  (i.e.  the  attitudes  intersect),  their  difference  is  the  4x1  zero  vector. 


qi(wi)  -qa(w2) 


0    0    0    0 


T 


(C.5) 


Basic  quaternion  multiplication  rules  can  be  used  to  express  qi(wi)  and  q2(^2)  in 
terms  of  the  components  of  the  vector  and  scalar  parts  of  the  true  attitude  quaternion. 

It  is  useful  to  consider  the  first  three  elements  of  the  attitude  ambiguity  quater- 
nions qi(wi)  and  q2(w2)  separately  from  the  fourth  element,  in  a  similar  manner  to  the 
partitioning  of  q.  These  first  three  elements  are  found  from  evaluation  of  equation  (C.2). 

qi(wi)i_3  =  cos(y)z  +  5sin(y)vi+sin(y)vi  x  z        i  as.  1,2  (C.6) 

where  x  represents  the  vector  cross  product.  The  fourth  elements  of  the  attitude  ambigu- 
ity quaternions  are  also  found  from  equation  (C.2). 

th(wt)4  =  scos(y)  -  sin(-rL)vi  •  z        i=  1,2  (C.7) 

where  (•)  represents  the  vector  dot  product. 

In  order  to  determine  the  intersection  of  these  two  attitude  ambiguities,  we  take  the 
difference  between  the  two  quaternions  to  find  which  wj  and  w2  cause  the  difference  to  be 
zero.  The  difference  of  the  first  three  elements  is  found  from  equation  (C.6). 


qi(^i)i-3  -q2(f>2)i-3  =  fcost—^z  +  ssint-^vj  +sin(y)v!  x  zj - 

cos(y )z  +  s  sin(y )v2  +  sin(y )v2  x  zj    (C.8) 
By  collecting  terms,  equation  (C.8)  can  be  written  as 

qi(^i)i-3  -  q2(^2)i-3  =  oz  +  svd  +  vd  x  z  (C.9) 
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where  we  have  defined  a  and  vd  as 


a      =      COs(y  )  -  COs(y  ) 

\d    =    sin(y)vi-sin(y)v2 


(CIO) 

(an) 


Similarly,  the  difference  of  the  fourth  components  of  qi(wi)  and  q2(^2)  may  be 
expressed  in  terms  of  the  new  parameters  a  and  vd. 


qi(wi)4-q2(w2)4    = 


wcos(y)-sm(y)vi  -v 


scos(y)z-sin(y)v2-v 


sin(y  )vi  -  sin(y  )v2 


+ 


=    — vj  •  v  +  as 


,wk         ,waxl 

COs(y  )  -  COs(y  )j 


(C.12) 
(C.13) 


The  attitude  ambiguity  intersections  are  found  where  (C.9)  and  (C.13)  are  equal  to 
zero.  The  free  parameters  in  these  two  equations  are  a  and  v^,  which  are  in  turn  set  by 
(Ji  and  u>2.  We  first  consider  the  situation  where  equation  (C.9)  is  set  to  the  3x1  zero 
vector.  Notice  that  of  the  three  vectors  in  equation  (C.9),  two  of  them,  namely  z  and  vd 
are  (by  definition,  always)  co-planer,  regardless  of  the  constants  that  scale  them.  The 
third  vector,  v  x  vd  is  always  orthogonal  to  the  plane  that  contains  the  first  two  vectors. 
So  for  equation  (C.9)  to  equal  zero,  the  sum  of  two  vectors  lying  in  a  plane  and  a  vector 
orthogonal  to  this  plane  must  equal  zero.  This  can  only  occur  if  the  two  coplanar  vectors 
are  equal  in  magnitude  and  anti-parallel,  and  the  length  of  the  orthogonal  vector  is  zero. 
For  this  particular  case,  the  length  of  the  orthogonal  vector  will  go  to  zero  when  the  two 
co-planar  vectors  are  parallel  or  antiparallel,  since  it  is  defined  as  the  cross  product  of  the 
coplanar  vectors.  So  for  a  solution  to  exist,  the  two  co-planar  vectors  must  be  equal  in 
magnitude  and  of  opposite  sign,  or  of  zero  length.  That  is,  either 


av  =  —  wvd 


(C.14) 
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or  both  a  and  vd  have  a  magnitude  of  zero.  The  consequences  of  the  condition  where  a 

and  vj  are  zero  is  discussed  later,  so  consider  the  case  where  equation  (C.14)  holds.  This 

may  be  written  as 

vd  =  — z  (C.15) 

s 

At  an  intersection  of  the  two  attitude  ambiguities  the  scalar  equation  (C.13)  will 
equal  zero  as  well.  Inserting  the  above  expression  for  vd  into  equation  (C.13)  yields 

0    =    -— z-z  + as  (C.16) 

s 


a 


1 

-  +  s 
s 


(C.17) 


which  has  a  solution  only  at  a  =  0.  Since  a  =  0,  equation  (C.14)  implies  that  vd  must 
be  zero  length  as  well,  since  s  is  arbitrary1    From  the  definition  of  a  and  \d,  and  the  fact 
that  Vx  7^  v2,  this  implies  that  at  an  intersection  of  the  two  attitude  ambiguities: 

uului  =  n  (2n)  (C.18) 

and  therefore  qi^)  and  q2(^)  are  equal  to  the  true  attitude  quaternion  q.  That  is, 

qi(wi)nq2(u;2)  =  q  (C.19) 

■ 
An  extension  of  this  to  multiple  satellites  is  straightforward.  Since  the  attitude  ambi- 
guities from  any  two  satellites  intersect  only  at  the  true  attitude,  the  intersection  of  all 
ambiguities  is  at  the  true  attitude  as  well. 


1  Even  if  s  =  0,  equation  (C.9)  can  only  equal  zero  when  a  =  0,  since  in  this  case  Izl  ^ 
0. 
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